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ABSTRACT 

We present a robust method to constrain average galaxy star formation rates, star formation histories, and 
the intracluster light as a function of halo mass. Our results are consistent with observed galaxy stellar mass 
functions, specific star formation rates, and cosmic star formation rates from z = to z = 8. We consider the 
effects of a wide range of uncertainties on our results, including those affecting stellar masses, star formation 
rates, and the halo mass function at the heart of our analysis. As they are relevant to our method, we also 
present new calibrations of the dark matter halo mass function, halo mass accretion histories, and halo-subhalo 
merger rates out to z = 8. We also provide new compilations of cosmic and specific star formation rates; more 
recent measurements are now consistent with the buildup of the cosmic stellar mass density at all redshifts. 
Implications of our work include: halos near 10 12 M Q are the most efficient at forming stars at all redshifts, 
the baryon conversion efficiency of massive halos drops markedly after z ~ 2.5 (consistent with theories of 
cold-mode accretion), the ICL for massive galaxies is expected to be significant out to at least z ~ 1-1.5, 
and dwarf galaxies at low redshifts have higher stellar mass to halo mass ratios than previous expectations 
and form later than in most theoretical models. Finally, we provide new fitting formulae for star formation 
histories that are more accurate than the standard declining tau model. Our approach places a wide variety 
of observations relating to the star formation history of galaxies into a self-consistent framework based on the 
modern understanding of structure formation in ACDM. 

Subject headings: dark matter — galaxies: abundances — galaxies: evolution — methods: N-body simulations 
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1. INTRODUCTION 

Constraining the buildup of stellar mass in galaxies pro- 
vides fundamental constraints on galaxy formation models. 
An ever growing number of galaxy observations have been 
made, cov ering a time perio d from 500 Myr after the Big 
Bang (e.g. IZheng et al.|[2012l) to the present day. A model 
of galaxy formation that can match observations over this en- 
tire stretch of time would represent a significant aid to our 
understanding. So far, matching the evolution of the stel- 
lar masses and star formation rates of galaxi es over this en- 
tire epoch with has proved difficult (see, e.g..lLu et al.M2012h 
iBoreani & Kravtsovl l20Tll: IWeinmann et all 120121) . Despite 
challenges of reproducing these observations with detailed 
models of galaxy formation, significant progress has been 
made in recent years with empirical models that connect the 
evolution of galaxy properties to the evolution of dark matter 
halos. 

Over the past decade, a range of studies have associated 
galaxies with dark matter halos at a given epoch, using a 
variety of techniqu es, including Halo Occupation Distribu- 
tion m odeling (e.g.. lBerlind & Weinberg||2002l : iBullock et alj 
20021) the Condi tional Luminosity Function modeling (e.g., 



Yang et aT] I2003[). and variants of the abundance matching 
technique (e! g ICohn et all \l99<k IKravtsoy & Klypinl [l999t 



Nevrinck et all 12004 IKravtsoy et all 12004 I Vale & Ostrikerl 



2004 I2006t Tasitsiomi et alj 12004 IConrov et al J 1 2006 1 : 
Shankar et all 120061: iBerrier et all 120061 iMarfn et alj booa 
Guo et alj 120 lOt IMoster et all 120 lOt iBehroozi et al.l 120 101) . 

The simplest abundance matching models assign the most 
massive observed galaxies in rank order to the largest halos in 



an equal simulation volume. With only slight modifications 
(e.g., using the peak mass for satellite halos), this technique 
accurately reproduces the redshift- and scale- d ependent clus- 
tering of galaxies (e.g., IConrov et al.l 120061 : iReddick et alj 
I20L1 . 

Because the cosmological model provides a prediction for 
the buildup of dark matter halos, one can combine knowl- 
edge of the galaxy-halo connection at different epochs with 
knowledge of the mass accretion and merger histories of ha- 
los to constrain the buildup of stellar mass in galaxies over 
cosmic time. This was firs t done in a comprehensive way by 
IConrov & Wechsler! (120091) . who provided an empirical con- 
straint on the star formation histories and stellar mass growth 
in galaxies from z = 2 to the present. This approach has 
been e xplored in several studies using variants of the tech- 
niques (IZheng et al.ll2007bt iWhite et alj |2007l: lLeitoeill20T2l: 
I Wang et alJ2012t IMoster et alj|2012l) . which together provide 
important constraints on the basic picture for the buildup of 
stars in galaxies and their connection to dark matter halos. 

These studies represented important advances, albeit with 
shortcomings. For example, most of these studies only mod- 
eled galaxy evolution from z ~ 2 to the present due to per- 
ceived conflicts between the integrated cosmic star formation 
and the cosmic stellar mass den sity at z > 1 dWilkins et alj 
12001: iHook ins & Beacomll200"6l) . In addition, because pas- 
sive stellar mass loss (e.g., supernovae of massive stars) de- 
pends on galaxies' star formation histories, these studies had 
to make assumptions about the ages of stars already formed 
at z > 2. Finally, these studies paid limited attention to the 
vast array of observational uncertainties as well as modeling 
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uncertain ties affecting their derived c onst raints; indee d , mos t 
results in lConrov& Wechsleri d2009t) and lWang etafl ( 12012b 
are presented without error bars. 

Here we present a new technique and a new compilation 
of observations, aimed squarely at resolving these issues. On 



high redshifts (e.g. IBouwens et al.ll201 It 


Bradlev et al.ll2012t 


IBouwens et al.ll2012l IMcLure et al.ll201 1 


) have made it pos- 



sible to place constraints on galaxy formation all the way to 
z ~ 8. At lower redshifts, we present the first constraints from 
stellar mass functions based on the PRIMUS survey, as well 
as fr om a new analysis at z = based on SDSS and GALEX 
data dMoustakas et alj|2012b . These are obtained from a con- 
sistent methodology from z = to z = 1, and thus allevi- 
ate many concerns about matching inconsistently-identified 
galaxies at different epochs. Using these combined data sets, 
we find that observations of the evolution of galaxy star for- 
mat ion rates and stellar mas ses can now be re conciled (see 
also iBernardi et aill2010l and iMoster et ai]|2012b . 

Our new method constrains the galaxy-halo relation using 
observed galaxy stellar mass functions as well as specific star 
formation rates and cosmic star formation rates. As with pre- 
vious studies, we match observed galaxies to halos, but the 
additional information on star formation rates allows us to 
break degeneracies and directly constrain the buildup of the 
intracluster light, as well as the amount of stars that can trans- 
fer from satellites to the central galaxy during mergers. We 
account for a number of statistical and systematic effects, in- 
cluding uncertainties from stellar population synthesis mod- 
els, dust models, star formation history models, the faint-end 
slope of the stellar mass function, observational complete- 
ness, scatter between stellar mass and halo mass. 

Given a parametrization for the intrinsic stellar mass — halo 
mass relationship as well as for these observational effects, 
we use a Markov Chain Monte Carlo (MCMC) method to 
find the allowed posterior distribution. This is important not 
only for addressing concerns that our results may be biased by 
limited observational constraints at high redshifts, but it also 
gives us the power to determine which observations would 
best improve the resulting constraints on the relationship be- 
tween stellar mass, star formation, and halo mass. The direct 
results of this analysis are constraints on the distribution of 
galaxy stellar masses as a function of halo mass and redshift. 
From these, many other constraints relevant to galaxy forma- 
tion are derived, including the average star formation rate as a 
function of halo mass and redshift, the average star formation 
history in galaxies at a given epoch as a function of galaxy 
stellar mass or host halo mass, the instantaneous baryon con- 
version efficiency of galaxies as a function of mass and red- 
shift, the buildup of the intracluster light, and the evolution 
of the stellar mass to halo mass ratio for progenitors of to- 
day's galaxies, all including full uncertainties and covering a 
redshift range from z = to z = 8. 

We provide a broad overview of our methodology in §|2] in- 
cluding our parametrization of the stellar mass - halo mass 
relation and the relevant uncertainties, with additional details 
in Appendices lAllBl We discuss the observational data sets 
relevant to our method in $3] with special attention to new 
measurements of cosmic star formation rates. We discuss 
the simulations that we use in $H along with recalibrations 
of the halo mass function (AppendixlGb, halo mass accretion 
rates (Appendix|H|, and subhalo merger/disruption rates (Ap- 
pendixQ}. We present our main results in $5] with discussion 
in $6] and a summary of our conclusions in $21 



Throughout this work, we assume a iChabrieri J20031) mi- 
tial mass function, the ste llar population synthesis model of 
Bruzual & Chariot d2003[) . and the dust model employed in 
Blanton & Roweisl (12007b . We convert data sets from other 
papers to these models as necessary. We additionally as- 
sume a fiat, ACDM cosmology with parameters CIm = 0.27, 
Q\ = 0.73, h = 0.7, n s = 0.95, and <j 8 = 0.82. 

2. METHODOLOGY 
2.1. Overview 

Much previous work has gone into determining the relation 
betw een halo mass and s t ellar mass as a function of redshift 
(e.g. IMoster et all l2012t lYang et all l2012t iLeauthaud et aT 



2012} iBehroozi et alJ 120101: IMoster et al.H2010t iGuoetal 
201(1: IConrov & Wechsleri 120091: IWang & Jin d l201oF 



Zheng et alj l2007bt iMandelbaum et al l)2006t iGavazzi et alj 



20071: I Yang et all l2009at lHansen et alj l2009fc iLin & Moh 
2004). The redshift evolution of this relation is due to three 
contributing physical effects for stars in galaxies (new star 
formation, merging satellite galaxies, and stellar mass loss) 
and one effect for halos (continuing mass accretion). Merging 
galaxies and mass accretion are well-constrained by dark 
matter simulations, and stellar mass loss is well-constrained 
with the assumption of an initial mass function (IMF) for 
stars. The most uncertain of these effects, in relation to dark 
matter halos, are new star formation and stars ejected into the 
intracluster light (ICL; or, at lower halo mass, the stars in the 
stellar halo) With observational constraints on the redshift 
evolution of the stellar mass and the star formation rate as 
well as computed constraints on satellite mass loss, we may 
effectively constrain both of these factors. 

Schematically, the star formation rate for a given galaxy 
over a given timestep is constrained by: 

SFR ■ At = Expected Stellar Mass Now 

- Remaining Stellar Mass from Previous Timestep 

- (Stellar Mass from Mergers) 

x(l -Fraction Ejected to ICL). (1) 

In this way, we can combine a specific assignment of stel- 
lar mass to halos with the evolution of the halo mass function 
over time, to derive star formation rates for halo trajectories. 
Our approach is to flexibly parametrize the possibilities for the 
stellar mass — halo mass relation, M„(M/,,z), as well as the 
uncertainties affecting the remaining terms in Eq. Q] Then, 
using a Markov Chain Monte Carlo method, we can deter- 
mine the allowed parameter space forM st (M/,,z) by comparing 
the implied SFRs and stellar mass abundances to observations 
across a wide range of redshifts. 

Parametrizing M*(M/,,z) separately for individual halos is 
beyond the scope of this work (although it will be explored 
in future papers). Here, we parametrize the median value of 
for halos at a given mass and redshift (i.e., a M*(M/,,z)) 
as well as the scatter in stellar mass at fixed halo mass as a 
function of redshift. By taking a summation of Eq. Q] over 
halos in a specific mass and redshift bin, we can derive the 
total amount of stars formed in that bin. Thus, dividing by 
the number of halos in that bin gives the average star for- 
mation rate, which is similar to approaches taken in previous 
work s (IConrov & Wechsl ell2009t iLel tnel 120121: IWang et alj 
l2012tlMoster etalj2012b . This approach is accurate except in 
terms of small second-order effects that come from different 
stellar mass histories having a range of stellar mass loss rates. 
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This aspect is discussed in detail in Appendix lAl Further dis- 
cussion of how individual galaxy star f orma tion histories can 
differ from the average is presented in 35.81 

In this section, we present our m ethod for parametrizing the 
stellar mass - halo mass relation ( 32.2. It , our para metriz ation 
of the uncertainties affecting obser vational data ( 32.2.21 ). our 
model for stellar mass los s (32.3K and our model for ejec- 
tion of stars into the ICL ( 32.41 ). To streamline the presenta- 
tion of this section, we present the methodology for incorpo- 
rating the effects of stellar mass accreted in mergers (which 
involves calibrating merger rates from simulations to z = 8) 
and for calculating observables (such as the stellar mass func- 
tion and specific star formation rates) to AppendiceslBlandlCl 
We present a summary of the methodol ogy and details of the 
Markov Chain Monte Carlo method in 



2.2. Determining the Stellar Mass - Halo Mass Relation to 

z = 8 

2.2.1. Intrinsic Relation 



As in iBehroozi et al.l ( 120101) . we parametrize the stel- 
lar mass - halo mass (SMHM) relation at a given epoch. 
The most commonly-used function i s a double power law 
(lYang et al.l 120 12b iMoster et all 120 121) . which has a charac- 
teristic halo mass (Mi) and stellar mass (M„,o), a low-mass 
slope (a), and a high-mass slope (/3). Because we include data 
down to low stellar masses (AT* - 10 7 ' 25 M Q ) and have tight 
observational constraints from the SDSS, we find that a dou- 
ble power-law cannot accurately fit the unique shape of the 
stellar mass function (SMF). As shown in AppendixlDl a dou- 
ble power-law results in a stellar mass function off by as much 
as 0. 1 dex at z = 0. For this reason, we choose a different, five- 
parameter form as the fitting function for the SMHM relation. 
The chosen form retains the low-mass power-law slope (a) as 
well as a characteristic stellar (M» o) and halo mass (Mi). The 
high-mass behavior is trickier to fit. We find the best match 
with a subpower (superlogarithm) function with index 7: 



Mi 

log 10 (M*(M* -> 00)) oc log 10 [ — 



(2) 



If 7 is 1, this equation becomes an ordinary power law. As 
7 approaches 0, this equation approaches a logarithm. In be- 
tween, this equation grows more slowly than any power law, 
but faster than any logarithm. While this equation can be- 
come multiply-valued for M/, < M\ , it is straightforward to 
mitigate this behavior as well as allowing a smooth connec- 
tion to a power-law forM/, < M\ . The specific parametrization 
we adopt is[] 

log 10 (M,(M /l )) = log 10 (eM 1 ) + /^log 10 Q| j\ -f(0) 

This is a power law with slope a for Mi, <C M\ and a subpower 
law with index 7 forM/j ^ M\ . The characteristic stellar mass 
to halo mass ratio is e at the characteristic halo mass M\ . The 
maximum errors of this fit to the SMHM relation are about 
0.025 dex, or roughly four times better than a double power- 
law fit, as discussed in AppendixlDl 

This relation gives the median stellar mass M» for halos of mass M/,. 
Because of scatter in the SMHM relation, the inverse of this function does 
not give the average halo mass for a given stellar mass. 



We use the virial mass (as defined in iBrvan & Normanl 
119981) to define the halo mass of central galaxies, and for satel- 
lites, we use the peak pro genitor virial mass (M nPa k). As this 
paper was being prepared, R eddick et al.l (120 121) found that the 
peak v max for halos is an even better proxy for stellar mass, 
at least at z = 0.05, in agreement with physical expectations 
that the depth of the halo potential well before it is impacted 
by stripping should be most correlated with the galaxy mass. 
The use of M pea k is at present easier both conceptually (to 
consider mass ratios between galaxies and halos) and opera- 
tionally (because it allows the use of several previously cal- 
ibrated relationships in the literature), so we stick with this 
choice in the present work for simplicity. The use of M pea k 
may result in slight underestimates of the clustering an d the 
satellite fraction (see discussion in Re ddick et al.ll2012l for a 
comparison of these quantities in the local universe). How- 
ever, 35.61 demonstrates that the differences between our re- 
sults and the best-fit ting stellar mass-halo mass relation of 
lReddicketai1(l20T2h are small, and we expect that they are 
within the current systematic errors. 

The redshift scaling of the relation poses a unique chal- 
lenge. At z = 0, the parameter fits must be flexible enough 
to match the tight constraints from SDSS observations. At 
higher redshifts, the fit must be flexible enough to allow 
exploration of possible star formation histories allowed by 
higher uncertainties, but not so flexible that it allows signif- 
icant over-fitting of the systematic and statistical errors in the 
data. Moreover, the scaling must be physical: obviously, neg- 
ative star formation rates cannot be allowed, nor can stellar 
mass to halo mass ratios that exceed the halo's baryon frac- 
tion. These physical constraints apply even at redshifts z > 8 
where we do not have data. 

By necessity, the error bars we obtain will be sensitive to 
the choice of the redshift fit. This is not only so for the varia- 
tion allowed in the stellar mass histories, but it is also true for 
the variation allowed in the nuisance parameters for system- 
atic errors. Because the latter are often poorly constrained by 
available data, the flexibility of the redshift fit also influences 
how much of the available nuisance parameter space can be 
explored. For this reason, we have tried a large number of 
different redshift scalings, including cubic spline interpola- 
tion with many control points (which results in over-fitting) 
and many different choices of scaling parameters with a or 
z — many of which have unphysical behavior at very high red- 
shifts. We have also tried fitting directly to abundance match- 
ing results (Appendix |E), but no robust trends with redshift 
emerged. 

The choice that we settled on has several conventions. For 
many variables in Eq.|3] we have three parameters: one each 
for the low-redshift value, one for the scaling at intermedi- 
ate redshifts (0.5 < z < 2), and one for scaling at high red- 
shifts (z > 2). We also apply an exponential shutoff to the red- 
shift scaling towards lower redshift to isolate the highly con- 
strained low-redshift parameters from higher-redshift ones. 
Our final parametrization is the following: 

v{a) = exp(-4a 2 ) 
log 10 (Mi ) = Mi.o + (Mi , fl (a - 1 ) +Mi r z)v 
log 10 (e)= e + (e a (a- l) + e : z)v + e a . 2 (a- 1) 
a = ao + (a a (a- \))v 
S = S Q + (S a (a-l) + d z z)v 
7 = 7 + (7 fl (a-l) + 7,7) l ,. (4) 
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Additionally, there will be scatter in stellar mass at fixed 
halo mass (£). While current studies have indicated that this 
is approximately 0.16-0.2 dex at z = ( iMore et al.l 120091: 
lYang et all l2009bt iReddick et all 120121) with no evidence 
for a mass trend at least down to 1O 12 M0 in halo mass 
(IReddick et al.ll2012l) . little is known about the redshift evolu- 
tion of this scatter. For this work, we parametrize the possible 
evolution of scatter with redshift via a two-parameter scaling: 



£ = 1), 



(5) 



and we ta ke the prior on £p to be 0.20 dex ±0.03 dex, as deter- 
mined by [Reddick et al.l (120121) using a combination of con- 
straints from the correlation function and conditional stellar 
mass function. 

2.2.2. Observational Systematics 

A full discussion of the systematic uncertainties affectin g 
stellar mass functions may be found in lBehroozi et al.l (120101) . 
These include uncertainties from the Initial Mass Function, 
the stellar population synthesis model, the dust model, the star 
formation history model, sample variance, Eddington bias, 
redshift errors, and magnification bias. Most of these effects 
result in a constant systematic offset (/x) in stellar masses: 



lQ gio 



M, 



A/ 



: /'- 



(6) 



However, some effects can have different behavior as a func- 
tion of stellar mass, largely on account of assumptions for 
passive galaxies (e.g., declining star formation rates and low 
dust fractions) not being true for active galaxies. We account 
for these effects by introducing an additional parameter (k) to 
capture offsets in the stellar masses of active galaxies: 



/ ^^*,meas,active . 

lo §i0 T7 I = jl+ K. 



M, 



*, true, active 



(7) 



For the fraction of galaxies that are qu iescent, we use a fitting 
formula for recent measurements from lBrammer et all d201 ll) . 
corrected to the stellar mass estimates we use in this paper: 



,/*passive(^*,meas j Z) ~ 



10 10.2+0.5 ZM 



-1.3 



+ 1 



(8) 



At low redshifts, most galaxies below 10 Mq in stellar 
mass are active; beyond z = 3, nearly all galaxies are consid- 
ered active. 

As in iBehroozi et al.l (120101) . we do not model uncertain- 
ties in the IMF in this work. The most significant remaining 
uncertainties come from the SPS model (~ O.ldex), the dust 
model (~ O.ldex), and the star formation history (~ 0.2dex). 
By definition, passive galaxies have had no recent star forma- 
tion history, so that the star formation history uncertainties are 
substantially smaller. For that reason, we take the priors on [i 
to be a log-normal distribution centered at zero with width 
0.14 dex (combined SPS and dust model errors), but the pri- 
ors on k to be of width 0.24 dex (all three sources combined). 
The functional form of how fi and k evolve with redshift is 
unknown. We therefore adopt the following fiducial formu- 
lae: 



^=^o+(fl-l)Ma 



(9) 
(10) 



At higher redshifts (z — 2 ,3). lMuzzin etalj ( 120091) suggest that 
uncertainties from SPS modeling become more important, on 
the order of 0.2 dex instead of 0. 1 dex. We thus set the widths 
of the priors on \x a and K a to be 0.22 dex and 0.30 dex, respec- 
tively, which take these increased uncertainties into account. 

Because stellar populations cannot be fully constrained 
with limited photometric information, stellar mass estimates 
for individual galaxies have intrinsic scatter relative to the 
true galaxy stellar mass . This causes an Eddington bias 
(lEddingtonl I1913L 119401) in the stellar mass function: on 
the high-mass end of the SMF, many low-mass galaxies are 
upscattered, but there are a limited number of higher-mass 
galaxies that can be downscattered. This results in a net in- 
crease in the observed numbers of high stellar-mass galax- 
ies. This effect is best estimated for each redshift range and 
modeling technique individually; while few authors have cor- 
rected for this error, many at l east p rovide an estimate of its 
effect. As in IBehroozi et al.l (120101) . we find that these es- 
timates depend most significantly on the redshift. We thus 
model the distribution in the observed stellar mass estimates 
(compared to the true stellar masses) as a log-normal Gaus- 
sian with mean and a redshift-dependent standard deviation 
given by: 

a(z) = a + a z z. (11) 



IConrov et alj 
= 0.07 and 



Following 
we take <tq 
a- = 0.04 ± 0.015, 
literature (IConrov et a 



Perez-Gonzalez et al 
Stark et all 120091: 



Marchesini et al.ll2Q10t) 



d2009b : IBehroozi et aT] J20T0I) . 
we take the prior on a- to be 
con sis tent wit h estimates fr o m the 
12001 iKaiisawa et al] [2009 : 
iMarchesini et all 2009 : 
iLee et all boi ll 



2008: 
Caputi et al.l 1201 It 



Finally, there are some concerns that, due to the Lyman 
break techniques used to detect high-redshift galaxies, not 
all high-redshift galaxies may be detected. If galaxies have 
bursty star formation, or if some fraction of star forming 
galaxies have extreme quantities of dust, then not all galaxies 
will be detected in Ly man break surveys (see also lStark et al] 
120091: Eee et alj|2~009l) . We parametrize the stellar mass com- 
pleteness of high-redshift surveys with a two-parameter fit, 
governing the amplitude (A) and redshift onset (z c ) of incom- 
pleteness: 

Cl(X 1 exp(z c -z)+l' 

Because there is little reason to believe that surveys are miss- 
ing a significant fraction of galaxies at z < 1, we require 
Z c > 0.8, and we set the completeness fraction to 1 for z < 1. 
This amounts to adopting the following formula for the com- 
pleteness: 



c(z) 



1 if z < 1 

C/ (z) + (l-c,(l)) if z > 1 



(13) 



We parametrize the effect of this incompleteness on the ob- 
served SFRs by a single parameter b, which sets the fraction 
of "bursty" vs. "dusty" star formation. Fully bursty star for- 
mation {b = 1) would have no effect on the total observed cos- 
mic SFR, but the observed SSFRs would be boosted by the 
incompleteness factor compared to the SSFR averaged over 
longer periods of time. Fully dusty star formation (b = 0) 
would lower the observed cosmic SFR by the same factor 
as the cosmic SM density. Presumably, a moderate amount 
of unobscured star formation would be enough to render a 
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high-redshift galaxy observable — hence, dustiness will im- 
pact galaxy completeness only if star formation is largely ob- 
scured. Given that models for SFRs from observed galaxies 
already attempt to correct for dust obscuration, the main effect 
of dustiness on SSFRs is that the SSFRs for very dusty galax- 
ies will not be observed. Therefore, in our current model, we 
assume that this has no effect on the average SSFR. 

2.3. Star Formation Histories 

Under our assumption of a IChabrierl {2003) IMF and 
iBruzual & Chariot! d2003i) stellar evoluti on tracks, we take the 
same model for stellar mass loss from IConrov & Wechslerl 
d2009h : 

(14) 



/ loss (f) = 0.051n 1 + 



0.3 Myr 

where /i oss (0 is the fraction of mass lost from a single stellar 
population as a function of the time since its formation. 

2.4. Intracluster Light and Merging Galaxies 

From Eq.Q] galaxies can build up stars either through merg- 
ers or internal star formation. Schematically, this implies that 
the unknown amount of stars deposited in mergers can be de- 
rived given separate constraints on the stellar mass growth of 
galaxies (i.e., from stellar mass functions) and their star for- 
mation rates: 



ASMdeposited = A5M centra i — SFR cenlr:1 \ At , 



(15) 



where ASMdeposited is the amount of stellar mass deposited 
from merging satellites, A5M cen trai is the stellar mass growth 
of the central galaxy, and SFR centm \ is the star formation rate 
of the central galaxy over the chosen period At. (Note that the 
real equation is slightly more complicated due to stellar mass 
loss from passive evolution; see AppendixlBl. 

When a galaxy merger occurs, the stars associated with 
the satellite galaxy may either be deposited onto the central 
galaxy or be ejected into the intracluster light (ICL). Because 
we can constrain the total amount of stars in merging satel- 
lites through the halo merger rate (AppendixUi and the stellar 
mass - halo mass relation, we can also derive the amount of 
stellar mass deposited into the ICL (AICL): 



AICL = ASM. 



incoming 



■ASM. 



deposited ■> 



(16) 



where ASMi ncom i ng is the total amount of stars in satellites that 
merge within At. 

This simple picture is made more complicated by inconsis- 
tencies in how galaxies and the ICL are separated. The total 
luminosity recovered and the fraction of it attributed to the 
galaxy are both dependent on surface-brightness limits, sky 
subtraction methods, and galaxy fitting methods, which are 
different for the different data sets we use and for different 
redshift ranges in the data sets themselves. Over most mass 
ranges and redshifts, this is not a problem: the rate of star 
formation in many galaxies results in an overall change in 
stellar mass between redshifts that is much larger than that 
attributable to a change in the luminosity modeling method. 
However, for massive galaxies (M» > 10 n M Q ) at low red- 
shifts (z < 1), this is a more significant issue. For example, 
including 5% more light between z = 0.3 and z = 0.1 in a 
massive galaxy (0.05 magnitudes) would result in an equal 
buildup in stellar mass as would a specific star formation rate 
of about2.5 x 10 -11 yr -1 (the expected SSFR for such galaxies; 
ISalim et al.ll2007l) . As the true SSFR drops to 10" 12 yr 1 , the 



galaxy models must be self-consistent to about 0.2% (0.002 
magnitudes) to avoid introducing a comparab le error, which 
is we ll beyond current calibration methods dBernardi et al.l 
l20T(ih . 

This said, changes in the mass of the galaxy because of defi- 
nitional issues are completely degenerate with incoming mass 
from mergers in Eqs. fl31and[T6T ASMdeposited is replaced by 
the sum of the definitional and deposited stellar mass changes 
(A5M d efinitionai+ ASMdeposited) everywhere it appears. Notably, 
Eq. [16] still gives us a way to robustly determine the buildup 
of stars in the ICL. 

If we used Eqs. [TBI and [TBI directly, the star formation rate 
would have to be parametrized separately from the stellar 
mass growth. Instead, we rewrite Eq.[l5]in terms of the frac- 
tion of stellar mass growth that comes from the star formation 



rate: 



ASMdeposited + ASMdefinitional = ( 1 ~ fsFlt) ASM C entral (17) 

S/^centralAf 



fsFR - ■ 



ASMcential 



(18) 



This is useful because there are strong physical priors on what 
fsFR can be. Low-mass galaxies do not have enough incom- 
ing stellar mass in satellites to account for a significant part 
of their stellar mass growth, so f SFR asymptotes to 100% for 
such galaxies. On the other hand, the highest- mass galax- 
ies ex perience almost no internal star formation ( S alim et al.l 
l2007h . meaning that /sfr must approach for such galaxies. 
So, we can approximate the halo mass dependence of f SFR 
using a double power law as 



fsFR(M h ) ■■ 



ICL 



+ 1 



(19) 



We allow Mj x i CL to be redshift-dependent using a two- 
parameter fit: 



M hJ cL = MucLfi + (a~ l)M hJCL . a . 



(20) 



Rather than add extra parameters for (3, we make the assump 
tion that /sfr(1O i6 M0) is 1%. This is equivalent to the con 
straint 

log 10 (99) 



16-log I0 (M/,. /C L)' 



(21) 



2.5. Methodology Summary 

Although the equations involved are somewhat compli- 
cated, the logical steps involved in our method are straightfor- 
ward and are shown visually in Fig.Q] These steps include: 

1 . Choose parameters for the stellar mass - halo mass re- 
lation as well as observational and methodology uncer- 
tainties (see Table [JJfor full list) using a Markov Chain 
Monte Carlo method. 

2. Use the chosen stellar mass - halo mass relation to pop- 
ulate halos with galaxies, and use merger rates and mass 
accretion histories calculated from dark matter simula- 
tions to infer galaxy growth rates. 

3. From the previous step, calculate the average star for- 
mation rates as a function of halo mass and redshift, 
and use those to calculate the average specific star for- 
mation rate and the cosmic star formation rate. In addi- 
tion, using the abundances of halos, calculate the stellar 
mass function. 
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l. Choose a stellar mass - 
halo mass (SMHM) relation 
from parameter space. 




2. Find galaxy growth histories O *(Z} 
by applying the SMHM relation < gb==-~0 
^ to dark matter merger trees. o '' o^"*^)' 




3. Derive the inferred stellar 
mass functions and star 3 
formation rates. 




Pi 




SSFR 






SM 


z 


SM 


4. Apply effects to simulate 
observational errors and & 
biases. 








SSFR 






SM 


z 


SM 



5. Compare to data and 
calculate likelihood of the '~ 
chosen SMHM relation. 




I 




SSFR 




SM 


z 


SM 



FIG. 1 . — Visual summary of the methodology used to constrain the stellar mass - halo mass relation. 

TABLE 1 
Table of Parameters 



Symbol 



Description 



Equation Section 



2.2.1 



M\ 
c 
a 
S 

7 
v 

t 

MhJCL 



Characteristic halo mass 
Characteristic stellar mass to halo mass ratio 
Faint-end slope of SMHM relation 
Strength of subpower law at massive end of SMHM relation 

Index of subpower law at massive end of SMHM relation 
Exponential cutoff of evolution of M„(M/,) with scale factor 
Scatter in dex of true stellar mass at fixed halo mass 
Characteristic halo mass at which half of stellar mass growth is due to mergers 



2.2.1 



2.2.1 



2.2.1 



2.2.1 



2.2.1 



2.2.1 

m 



Systematic offset in stellar masses and SFRs 
k Systematic offset in stellar masses for active galaxies 

(7 Scatter in measured stellar mass at fixed true stellar mass 

c Galaxy detection completeness 

b Fraction of incompleteness due to burstiness (as opposed to dustiness) 

p Correlation of SFR to stellar mass at fixed halo mass 



2.2.2 



2.2.2 



2.2.2 



2.2.2 



2.2.2 
"Cj 



4. Apply corrections for observational errors and biases 
to the derived stellar mass function and star formation 
rates. 

5. Compare with observational measures of the stellar 
mass function and star formation rates (i.e., sum x 2 er- 
rors for each data point) to calculate the likelihood that 
the chosen stellar mass - halo mass relation matches 
observed results (i.e., exp(-0.5x 2 ))- 

6. Return to step #1 until the MCMC algorithm has con- 
verged. 

To ensure convergence on a space with such a large num- 
ber of paramet ers, we use the A daptive Metropolis explo- 
ra tion method (fHaario et all 1200 ll) . Although the algorithm 
of iDunklev et all d2005h indicates that we converge after only 
5 x 10 5 points, we continue to run for 4 x 10 6 total points to 
minimize the chance of unexplored regions and to ensure that 
the adaptive updates of the step covariance matrix converge 
and do not bias our final results. 



3. OBSERVATIONAL CONSTRAINTS AND UNCERTAINTIES 

We use three types of observational constraints on star for- 
mation in galaxies: the stellar mass function (SMF), the spe- 
cific star formation rate of galaxies (SSFR, star formation 
rate per unit stellar mass), and the cosmic star fo rmat ion rate 
(CSFR). Stellar mass functions are dis cusse d in 33.11 cosmic 
star formation rates are discuss ed in 33.31 and specific star 
formation rates are discussed in 33.21 

To ensure a fair comparison between the different data sets, 
we convert each measurement to the same set of systematic 
assumptions. For the i nitial mass func tion (IMF), we con- 
vert all results to that of lChabrierl (120031) . We do not consider 
evolving IMFs, nor do we consider the uncertainties for a non- 
universal IMF, as both uncertainties are beyond the scope of 
this paper. For the stellar population synt hesis and dust mod- 
els, we co nvert all results to the mod els of iBruzual & Charloi 
(2003J and lBlanton & Roweisl (120071) . In all cases, we convert 
authors' supplied units to physical units under the assumption 
that h = 0.7. 
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TABLE 2 
Table of Priors 



Symbol 



Description 



Equation 



Prior 



£o Value of f at z = 0, in dex 

§ a Redshift scaling of £, in dex 

fig Value of /i at z = 0, in dex 

fi a Redshift scaling of fi, in dex 

kq Value of re at z - 0, in dex 

re a Redshift scaling of re, in dex 

<r, Redshift scaling of cr, in dex 

A Amplitude of galaxy detection completeness 

z c Onset redshift of galaxy detection completeness 
b Fraction of incompleteness due to burstiness (as opposed to dustiness) 

po.5 Correlation of SFR to stellar mass at fixed halo mass at a = 0.5 



G(0.16,0.04) 
G(0,0.16) 
G(0,0.14) 
G(0,0.22) 
G(0,0.24) 
G(0,0.3) 
G(0.05. 0.015) 
U(0, 1) 
Zc > 0.8 
U(0, 1) 
U (0.23, 1.0) 



TABLE 3 

Observational Constraints on the Stellar Mass Function 



Publication 



Baldrv et al. (2008) 
Moustakas et al. (2012) 
F^erez^ojvzalez_etaL Q00S^) 
Mortlock et al. f2011) 
Marchesini et al. (2009) 
Marchesini et al. (2010) 
Lee et al. (2012) 



Stark et al. (2009) 
Flojrwens_etaL Q0FT) 
Bradley et al. (2012) 



Redshitts 


Colors 


Area 


Notes 


0.003-0.05 


ugriz 


4783 deg 2 


VBD 


0.05-1 


UV-MIR 


9deg 2 


D 


0.2-1.6 


UV-MIR 


0.184 deg 2 


ISD 


1-3 


BVizHim 


0.0125 deg 2 


VID 


1.3-4 


B-MIR 


0.142 deg 2 


ID 


3-4 


UV-MIR 


0.43 deg 2 


ID 


4-5 


B-MIR 


0.089 deg 2 


SD 


6 


fi-MIR 


0.089 deg 2 


SD 


7-8 


B-Hi m 


0.0148 deg 2 


UD 


8 


MIR 


0.076 deg 2 


UD 



NOTE. — Letters correspond to the following corrections made to the 
published results: I (Initial Mass Function), S (Stellar Population Synthe- 
sis model), D (Dust model), V (Sample Variance), B (Surface Brightness 
incompleteness), C (Cosmology, specifica lly h), U (UV to stellar mass con- 
version accordin g to Gonzalez et al. 20TT|). The local results (z < 0.2) in the 
Moustakas et al. (2012) mass functions are taken from the Sloan Digital Sky 
Survey and cover an area of 2505 deg 2 . 



3.1. Stellar Mass Functions 



NOTE. — G(.v, y) denotes a Gaussian distribution with center x and width y. U{x\ ,X2) denotes a uniform distribution from x\ to #2. The remaining parameters 
have no explicit priors; M\ and e are explored in logarithmic space, whereas ce, 8, and 7 are explored in linear space. 

this redshift range come from several previously published re- 
sults that suggest that the evolution in the cosmic SFR density 
is inconsist ent with the estimated cosmic stellar mass den- 
: sitv at z > 1 (|Nagamine et al.l200ollHopkins & Beacoml2006l: 
- iPerez-Gonzalez et al.ll2Q08t [wilkins et al.ll2008l) . One possi- 
ble explanation is an evolving or a non-universal IMF; a num- 
ber of di fferent lines of evidence exist in support of this the- 
ory (e.g. .Ivan Dokkum & Conrovl2010tlLucatello et al.ll2005t 
iTumlinsonl 12007 allbt Ivan Dokkuml l200cC Another possible 
reason is that star formation at high redshift could be bursty 
or dust-obscured, yield ing an incomplete census of galax ies in 
Lvman-break s urveys dLee et al.ll2009HStark et al.ll2009l) . On 
the other hand. lReddy & Steidell (120091) offer a simpler expla- 
■ nation. They appeal to luminosity-dependent reddening cor- 
rections in the ultraviolet luminosity functions at high redshift 
as well as steeper faint-end slopes for stellar mass functions, 
and demonstrate that the purported discrepancy then largely 
vanishes. 

In our analysis of recent cosmic SFR literature ( 93.2t . we 
find that newer estimates of the cosmic SFR are systemat- 
ically lower at z > 3 than the data in iHopkins & Beacoml 
d2006l) (see F i g. |2). A similar conc l usion is reached in 
iBernardi et"ail d2010l) and iMoster et al.l ( 12012b . In addition, 
deep probes of the stellar mass function at lower redshifts 
(z < 1 -5) hav e resulted in findi ngs of steep faint-end slope s 
(Baldry et al.ll2008b iDrory et alJEOOfl iMortlock et aDl201 lh . 
We fin d this to be strong eviden ce in support of the explana- 
tion in iReddy & Steidell d2009). Nonetheless, we also con- 
sider simple models for the effects of bursty or dusty modes 
of star formation; an analysis of models for an evolving or 
non-universal IMF is beyond the scope of this paper. 

The highest redshifts present unique challenges. Herein, 
for redshifts 7 < z < 8.5, we convert the UV luminosity func- 
tion given in iBouwens et al.l (120111) to a stellar mas s func- 
tion according to the recipe in iGonzalez et al.l d201 lh . It has 
been suggested that this conversion is reasonable given the 
blue UV continuum slopes of such high-redshift galaxies, 
and follow-up observations from I RAC that imply th at these 
early galaxies are nearly dust-free (lLabbe et al.ll2010h . At the 
same time, conc erns arise because the UV luminosity func- 
tion presented in IBouwens et al.l (1201 ll) is steep (a pa -1 .6 to 
-2.0, instead of the shallower a w — 1.5 to -1.7 at lower red- 
shifts). One interpretation of this result is that assumptions 
for galaxy morphology (and, therefore, the ability to mea- 
sure estimate the full UV flux or half-light radii) can bias 
the calcul ated UV luminosity f unction to have a steeper faint- 
end slope dGrazian et al.ll20lH) . Indeed, any systematic effect 
that interferes more with brighter galaxies than with dimmer 



We combine many overlapping data sets dBaldrv et al.l 
20081: iMoustakas et al.1 120121: IPerez-Gonzalez et al.l 120081: 
Mortlock et all 120111: iMarchesini et al.l I2009t iStarketafl 
20091; IBouwens et al.ll201 ll) to constrain the evolution of the 
stellar mass function from z = to z = 8. These are shown in 
Table |3] 

Different observational challenges present themselves in 
each redshift range. At low redshifts, the statistics are suf- 
ficient that systematic errors in ca lculating stellar masse s are 
the largest source of uncertainty dBehroozi et al.l l20Toh . At 
higher redshifts, this is no longer the case; statistical and sam- 
ple variance errors can be of equal magnitude to estimated 
systematic errors. Nonetheless, under the assumption of a 
smoothly varying stellar mass - halo mass (SMHM) relation, 
our approach (which simultaneously constrains the SMHM 
relation against many stellar mass functions at different red- 
shifts) can constrain the statistical errors better than the indi- 
vidual error bars on stellar m ass functions m i ght su ggest. We 
refer the interested reader to lBehroozi e t al. (2010]) for a dis- 
cussion of the most common uncertainties affecting observa- 
tion constraints on the stellar mass function and to §12.2.1l for 
our parametrization of their effects. In this section, therefore, 
we limit ourselves to discussing special concerns relevant to 
using stellar mass functions over this broad range of redshifts. 

Questions about the accuracy of stellar mass estimates over 
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ones (e.g., dust) will tend to cause a n overestimation of the 
slope of the UV luminosity function. iGonzalez et al.l d201 lb 
assume that luminosity-dependent corrections for the UV-to- 
stellar-mass conversion do not evolve substantially from z = 4 
to z = 7. This may in fact overestimate the dust corrections, 
an error that would lower the faint-end slope, partially cor- 
recting for a failure to account for galaxy morphology. As the 
use of UV-only stellar masses therefore remains controver- 
sial, we presen t resul ts both with and without the converted 
iBouwens et al.l (1201 ll) stellar mass functions for z = 7 - 8 in 
Appendix|J] 

3.2. Cosmic Star Formation Rates 

To assemble constraints on cosmic star formation rates, we 
conducted a comprehensive literature search on the astro-ph 
arXiv in September 201 1; we selected papers that 

1. Were posted within the past 5 years (2006-201 1). 

2. Contained the word "cosmic" and either the word 
"SFR" or the words "star" and "formation" in the title 
or abstract. 

3. Contained an estimate of the cosmic star formation rate 
from observations. 

4. Were not superseded by later publications. 

5. Contained a reference to the initial mass function (IMF) 
used. 

The selected papers are summarized in Table @] Many dif- 
ferent types of surveys were conducted, including estimates 
of the cosmic SFR from narrowband (e.g., Ha), broadband 
(UV-IR), and radio (1.4 GHz) surveys. 

As shown in Fig. [2 newer estimates of the cosmic SFR are 
systematically and substantially lower tha n previous compi- 
lations (such as iHopkins & Beacomll2006l) at z > 3. At high 
redshifts, the predominant method of probing star formation is 
through U V emission; however, conversion to a star formation 
rate re quires quantifying dust obs curation systematics. Previ- 
ously, IHopkins & Beacoml d2006l) assumed fixed dust obscu- 
ration for z > 2.5. However, new measurements have shown 
that the amount of dus t is likely to decrease substantially with 
redsh ift beyond z ~ 3 (IBouwens et al.ll20 1 2h iReddy & Steidel 
120091) . The inferred star for mation rates thus were like ly to 
have been overcorrected in IHopkins & Beacoml ( I2006T) . al- 
though substantial uncertainty about the cosmic star forma- 
tion rate at high redshifts still remains. 

Because most survey authors do not fully consider sys- 
tematic errors, most reported error estimates (1-10%) are too 
small to explain the observed variance in published estimates 
of the cosmic star formation rate (see Fig. [2]). Therefore, in- 
stead of using authors' estimates for systematic uncertainties, 
we estimate the true systematic errors by computing the vari- 
ance of our collection of published SFR estimates, as detailed 
in Appendix [F] Our results are given in Table [5] We find the 
average systematic error to range from 0.13 dex at z = to 
0.27 dex for z> 3. 

3.3. Specific Star Formation Rates 

We conducted a literature survey in an almost identical 
manner as for the cosmic star formation rate to assemble a 
collection of specific star formation rate (SSFR) constraints 
as a function of stellar mass, again in September 2011. In 
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FIG. 2. — Cosmic star formation rates reported in the past five years have 
been dramatically lower at high redshifts (z > 2) than previous results. Top 
panel: our new compilation of cosmic star formation rates, which excludes 
those reported before 2006 (see Table [4}, along with a best-fitting double 
power la w. Bottom panel: the compi lation of cosmic star formation rates re- 
ported in ( Hopkins & Beacom 2006) along with a best-fitting double power 
law (black dashed line) compared to the best-fitting double power law from 
our new compilation (blue dashed line). The largest differences at high red- 
shifts result fro m c hanging assumptions about the amount of dust present at 
z > 3. Fits (Eq.|FTJ to the new and older data sets are given in Table |7J Error 
bars are not shown since publications have often dram atically underestimated 
the true magnitude of the systematic errors (see j)3.2) . 



this case, the selection criteria on astro-ph were papers that 
were posted within the past 5 years (2006-2011); contained 
the word "specific" and either the word "SFR" or the words 
"star" and "formation" in the title or abstract; contained an 
estimate of the specific star formation rate as a function of 
stellar mass, both as calculated from observations; were not 
superseded by later publications; contained a reference to the 
initial mass function (IMF) used; and did not limit themselves 
to a specific class of galaxies (sub-mm, ULIRGs, satellites, 
star-forming only, etc.). The selected papers are summarized 
in Table [6] 

Most authors do not consider the full range of systematics 
that could affect their stellar masses or SFR est imat es. To 
accout for this we take a similar approach as in 33.21 and es- 
timate systematic errors by calculating inter-publication vari- 
ances. Details of our approach are given in Appendix|F] We 
did not find evidence for mass or redshift trends in the vari- 
ance, and have adopted a uniform systematic uncertainty of 
0.28 dex for SSFRs. 

Besides these observational constraints, we apply one mod- 
eling constraint to SFRs in high-mass halos at low redshifts. 
The slope of the stellar mass — halo mass relationship flattens 
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TABLE 4 

Observational Constraints on the Cosmic Star Formation 
Rate 



TABLE 6 

Observational Constraints on the Specific Star Formation 
Rate 



Publication Redshifts Technique 



Redshifts 


Technique 


Area 


0.005-0.2 


UV 


741 deg 2 


0.2-1 


UV/IR 


0.458 deg 2 


0.2-1.1 


UV/IR 


0.5 deg 2 


2.2 


Ha 


0.0156 deg 2 


1.4-2.5 


UV- 1.4GHz 


<0.025 deg 2 


0.2-3 


1.4 GHz 


2 deg 2 


0.5-3.5 


UV/IR 


0.029 deg 2 


3.3-4.3 


UV/IR 


5.3 deg 2 


0.4-5 


UV/IR 


0.025 deg 2 


4-8 


UV 


0.040 deg 2 


6-8 


UV 


2 deg 2 


6-8.7 


UV 


0.0125 deg 2 



Area 



Notes 



Notes 



Robotham& Driver (2011) 0.013-0.1 UV 

Salim et al. (2007) 0.005-0.2 UV 

Lvetal. (2011a) 0.8 Ha 

Zheng et al. (2007a) 0.2-1 UV/IR 

Ruiopakarn et al. (2010) 0-1.2 FIR 

Smolcic et al. (2009) 0.2-1.3 1.4GHz 

Shim et al. (2009) 0.7-1.9 Ha 

Tadakietal. (2011) 2.2 Ha 

Haves et al. (2010) 2.2 Ha 

Magnelli et al. (2011) 1.3-2.3 IR 

Karimetal. (2011) 0.2-3 1.4 GHz 

Lvetal. (2011b) 1-3 UV 

Kaiisawa et al. (2010) 0.5-3.5 UV/IR 

Dunne et al. (2009) 0-4 1.4 GHz 

Cucciati et al. (2012) 0-5 UV 

Le Borgne et al. (2009) 0-5 IR-mm 

van der Burg et al. (2010) 3-5 UV 

Yoshida et al. (2006) 4-5 UV 

4-8 UV 



Bouwens et al. (2012) 



833 deg 2 
741 deg 2 
0.8 deg 2 
0.458 deg 2 
0.389 - 9 deg 2 
2 deg 2 
0.029 deg 2 
0.0156 deg 2 
0.0156 deg 2 
0.0786 deg 2 

2 deg 2 
0.241 deg 2 
0.029 deg 2 

0.8 deg 2 
0.611 deg 2 
varies 
4 deg 2 
0.243 deg 2 
0.040 deg 2 



NOTE. — Letters correspond to the following corrections made to the pub- 
lished results: I (Initial Mass Function), A (Intrinsic Scatter Correction; see 
AppendixIPY The technique of Le Borgne et al. 1 2009) (parametric derivation 
of the cosmic SFH from counts of IR-sub mm sources) uses multiple surveys 
with different areas. 



TABLE 5 

Error Estimates for the Cosmic Star Formation Rate 



Redshift Range 


Inter-Publication 


This Work - HBOft 


0.025-0.5 


0.13 


-0.08 


0.5-0.9 


0.13 


-0.12 


0.9-1.5 


0.17 


-0.07 


1.7-3 


0.19 


-0.05 


3-8.0 


0.27 


-0.33 



Salim etal. (2007) 
Zheng et al. (2007a) 
Noeske et al. (2007) 
Tadaki etal. (2011) 
Daddi et al. (2007) 
Karim etal. (2011) 
Kaiisawa et al. (2010) 
Lee etal. (2011) 
Feulner et al. (2008) 
Bouwens et al. (2012) 
Schaerer & de Barros (2010) 
Al McLure et al. (2011) 

Al 

NOTE. — The le tter "I" corresponds to an Initial Mass Function correction. 
IDaddi et alj 0001) report using a fractio n of the GOODS-N and GOO DS-S 
fields, but did not report an exact area. Schaerer & de Barros 12010) use a 
combination of fields, the largest of which is COSMOS (with 2 deg 2 area). 
However, this has a comparatively bright limiting magnitude. 



NOTE. — "Inter-publication" refers to the standard deviation in dex be- 
tween published data points (from references in Table f4j in redshift bi ns af - 
ter correction for the evolution of the cosmic star formation rate (Eq. [FT}; 
the error in the error estimate for each redshift bin is of order 25%. "T his 
Work - HB06" refers to the difference in dex between t he best fit (Eq . lFTl to 
data p oints in Table [4] and the best fit to data points in Hopkins & Beacom 
12009) over the same redshift bins; this is a qualitative estimate of how much 
assumptions about calculating the cosmic star formation rate have changed 
over the past decade. 



towar ds high halo masses dBehroozi et al. 2 0~10t IMoster et all 
l2010h . meaning that central galaxy stellar mass becomes a 
poor indicator of the host halo mass. Consequently, ob- 
served SSFRs (reported as functions of stellar mass) offer 
poor constraints on the star formation rate of high-mass halos. 
Nonetheless, observations of galaxy clusters have indicated 
that the centr al galaxies generally have very low star forma- 
tion rates (e.g. lDonahue et al.l2010tlRawle et al.l2012l). To ac- 
count for this we apply a weak prior based on the lWetzel etail 
(12011) group catalog. The recent star formation history of 
Mh > 10 I4 M Q halos is averaged from z = to z = 0.2 (corre- 
sponding to the range of the SDSS survey used in lWetzel et al.l 
120 121) : if it is above 1M Q yr" 1 , the the total \ 2 error of the fit 
is increased by ((SFH) - 1M yr _1 ) 2 (M Q yr" 1 )" 2 . 

4. SIMULATION DATA 

We rely mainly on the high-resolution Bolshoi simulation, 
described in iKlvpin et all (1201 lb . Bolshoi follows 2048 3 



(?s 8.6 billion) particles in a comoving, periodic box with 
side length 250 h~ x Mpc from z = 80 to the present day. Its 
mass resolution (1.9 x 10 8 M Q ) and force resolution (1 h~ l 
kpc) make it ideal for studying the evolution of halos from 
10 10 M Q (e.g., satellites of the Milky Way) to the largest clus- 
ters in the universe (10 15 M Q ). Bolshoi was run as a colli- 
sionless dark matte r simulation with the Adaptive Refinement 
Tree Code (ART; iKravtsov eFaTI 1 1 997h IKravtsov & Klvpinl 
119991) assuming a flat, ACDM cosmology (f2 M = 0.27, f2 A = 
0.73, h = 0.7, erg = 0.82, and n s = 0.95). These cosmologi- 
cal parameters are con sistent with results from both WMAP5 
dKomatsuet al.l 120091) and the latest WMAP7+BAO+H re- 
sults dKomatsu et alJl201 ll) . 

Halos in Bolshoi were identified using ROCKSTAR, a seven- 
dimensional hal o finder that uses pha se space plus tempo- 
ral information dBehroozi et al.l 1201 lal) . Me rger trees were 
enerat ed using a new algorithm presented in iBehroozi et alj 
201 lbl) . which enforces physically consistent evolution of 
halo properties across timesteps to provide increased robust- 
ness against anomalies in the halo finder. Halo masses are cal- 
culated using spherical overdensities, according to the virial 
overdensity criterion of iBrvan & Norman] dl998l) . 

To supplement results from Bolshoi where there exist con- 
cerns about statistics (i.e., for high-mass halos) or for initial 
conditions (as Bolshoi uses Zel'dovich/ILPT initial condi- 
tions), we also make use of two larger simulations. MultiDark 
follows the same number of particles in a much larger volume 
(1000 h~ l Mpc) with 7 times worse force resolution and 64 
times worse mass resoluti on, but using the i dentical simula- 
tion code and cosmology dRiebe et alj 1201 ll) . Consuelo fol- 
lows 1400 3 particles in a somewhat larger volume (420 h 
Mpc) than Bolshoi; it has eight times worse force resolution 
and fourteen times worse mass resolution (McBride et al, in 
preparation; see also IBehroozi et al1l201 l b) PI It has a slightly 
different cosmology (fi,„ = 0.25, J1a = 0.75, h = 0.7, n, = 1.0, 
and a = 0.8) and was run using the GADGET-2 code dSpringell 
120051) . Notably, however, it uses the more accurate 2LPT ini- 
tial conditions instead of Zel'dovich initial conditions, so it 
serves as a way to correct for inaccuracies in the high-redshift 

2 |http : //lss .phy . vanderbilt . edu/ lasdamas/l 
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FIG. 3. — Top panel: Evolution of the stellar mass function from z = 
to z = 8 in the best fitting model (colored lines), compared to observations 
(points with error bars; for clarity not all data is shown). Bottom panel: 
Observational constraints on the cosmic star formation rate (black points), 
compared to the best-fit model (red solid line) and the posterior one-sigma 
distribution (red shaded region). 

mass function for Bolshoi. 

Results from these three simulations are used to calibrate 
the halo mass function (AppendixlGil. halo mass accretion his- 
tories (Appendix IHl. and subhalo disruption rates (Appendix 
IJl used in this work. 

5. RESULTS 

The method presented above results in a posterior distribu- 
tion for the set of parameters describing models that match 
observed stellar mass functions, specific star formation rates, 
and cosmic star formation rates from z = to z = 8. Our 
best-fitting parameters with one-sigma limits are as follows: 

Intrinsic Parameters: 

v = exp(-4a 2 ) 



X 10 



a- 

00 



V3 10 
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FIG. 4. — The best fitting model (red line) and posterior one-sigma distri- 
bution (red shaded region) for the evolution of the specific star formation rate 
from z = to z = 8, compared to observational estimates (black points). 
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log 10 (M VC L) = 12.399^|+ (-2.644^)(a-l) 



Our total x 2 error from all sources (observational and theoret- 
ical) is 235. For the number of observational data points we 
use (584), the nominal reduced \ 2 is 0.4. While the true num- 
ber of degrees of freedom is not easily computed (note for 
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FIG. 5. — Left panel: Average star formation rates as a function of halo mass and redshift. The overlaid white lines show average mass accretion histories 
for halos as a function of redshift for comparison. The grey area shows halos that would have a mass of > 10 15 5 Mq at z = and therefore are not expected to 
exist. Right panel: Star formation histories (SFH) as a function of present-day halo mass and redshift, for galaxies at z = 0. This figure shows the historical star 
formation rate for stars in the galaxy at the present day. Since the contribution of stars from merging galaxies is so low, this is equivalent to the star formation 
rate traced along the white mass accretion trajectories in the left panel. 
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FIG. 6. — Left panel: Average star formation rates for the galaxies in halos at a given halo mass and redshift (lines). Shaded regions indicate the one-sigma 
posterior distribution. Right panel: Average star formation histories as a function of halo mass and redshift (lines). Shaded regions indicate the one-sigma 
posterior distribution. Histories for 10 15 Mq halos are not shown as they are very similar to those for 10 14 Mq halos. 
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FIG. 7. — Left panel: Evolution of the derived stellar mass as a function of halo mass. In each case, the lines show the mean values for central galaxies. These 
relations also characterize the satellite galaxy population if the horizontal axis is interpreted as the halo mass at the time of accretion. Error bars include both 
systematic and statistical uncertainties, calculated for a fixed cosmological model (see j4]for details). Right panel: Evolution of the derived stellar mass fractions 
(M*/M;,) as a function of halo mass. 
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FIG. 8. — Evolution of the derived stellar mass fractions (M*(z)/M/,(z)) as 
a function of halo mass at the present day. More massive halos used to have 
a significantly larger fraction of mass in stars, but the peak star formation 
efficiency has remained relatively constant to the present day. 

example that covariance matrices are not available for most 
stellar mass functions in the literature), this suggests that our 
best fit is a reasonable match to the data. This is shown visu- 
ally in Fig. [3] which shows the evolution of the stellar mass 
function and the posterior distribution for the observed cosmic 
star formation rate and in Fig. @1 which shows the posterior 
distribution of specific star formation rates. 

As discussed in $2j each model in the distribution contains 
complete information on the average stellar content and for- 
mation history of halos as a function of mass and redshift. 
Herein, we focus on a few interesting implications, leav ing 
more extensive coverage to a follow-up paper. In 35.11 we 
present derived star formation rates and star formatio n hi sto- 
ries as a function of halo mass and redshift. Next, in 35.21 we 
discuss constraints on the stellar mass to halo mass relation, 
both as a function of h istor ical halo mass and as a function of 
halo mass at z = 0. In 35.31 we compare the trajectories of the 
stellar ma ss an d intracluster light (ICL) buildup in galaxies. 
Then, in 35.41 we discuss instantaneous baryon conversion 
efficiencies and how they relate to integrated baryon conver- 
sion efficiencies as well as stellar ages and formation times in 
35.51 In 35.61 we compare our main results to those obtained 
by previous studies. We discu ss th e main effe cts o f the uncer- 
tainties we have modeled in 35.71 finally, in 35.81 we present 
new fitting formulae for individual galaxy star formation his- 
tories relevant to observers. 

5.1. Star Formation Rates and Histories 

We show derived star formation rates as a function of halo 
mass in the left panel of Fig. [5] and the corresponding star 
formation histori es f or galaxies at z = in the right panel. 
As discussed in 35.31 the contribution from merging galaxies 
to the central galaxy is small, so star formation histories for 
galaxies at z = trace the star formation rate as a function of 
their progenitor's halo mass. Also, for this reason, star forma- 
tion histories for galaxies at a given redshift z g are nearly the 
same as for the the z = stellar populations, except truncated 
at z = z g U 

We show similar plots with one-sigma uncertainties in Fig. 

3 Note that this also requires that the halo accretion histories are similar 
between the z = z g progenitors of z = halos and all similar halos of the same 
m ass as the progenitors at z = z g ', however, this has been shown to be the case 
in lMcBride eiaU f2009h . 



[6] The left-hand panel demonstrates that the star formation 
rate at fixed halo mass has been monotonically decreasing 
since very early redshifts. This rate of decrease is different for 
different halo masses. At moderate to high redshifts (z > 2), 
larger halo masses generically have larger average star for- 
mation rates. However, at lower redshifts, the highest mass 
halos (Mi, > 10 14 M Q ) become so inefficient that they have 
lower star formation rates than group-scale (10 13 M Q ) halos 
or Milky-Way sized (10 12 M Q ) halos. 

From the perspective of individual galaxies, it is more il- 
luminating to look at the star formation history in the right 
panel of Fig. [6] Because halos continually gain mass over 
time, and do so more rapidly at early redshifts, the star for- 
mation history for galaxies is not monotonically decreasing. 
Instead, it increases rapidly with time, approximately as a 
power law in time. Depending on present-day halo mass, 
the galaxy's star formation rate reaches a peak at a redshift 
between z = 0.5 to z = 2.5 (higher redshifts for higher halo 
masses) and then decreases until the present day. The rate 
of decrease depends again on the halo mass, with high halo 
masses shutting off more rapidly than lower halo masses. 
Cluster-scale (M/, > 1O 14 M ) halos form most of their stars 
rapidly, at early times, whereas galaxies in Magellanic Cloud- 
scale halos (1O U M ) form stars over an extended period of 
time (see also 85.4) . 

5.2. The Stellar Mass - Halo Mass Relation 

We show constraints on stellar mass — halo mass (SMHM) 
relation from z = to z = 8 in the left panel of Fig. [JJ and on 
the stellar mass — hal o mass ratio in the rig ht panel. As seen 
in our previous work dBehroozi et al.ll2010l) . there is a strong 
peak in the stellar mass to halo mass ratio at around 10 12 M© to 
at least z ~ 4 and a weaker peak still visible to z ~ 8. While the 
location of the peak appears to mo ve to higher masses wit h in- 
creasing redshift (consistent with iLeauthaud et al. 2012[), the 
abundance of massive halos is also falling off with increasing 
redshift. 

In terms of dwarf galaxies (M/, ~ 10 10 M Q ), we only have 
constraints from observations at z = 0. These have been the 
subject of recent interest due to the finding of higher- than- 
expected stellar mass to halo mass ratios in dwarf galaxies 
around the Milky Way dBovlan-Kolchin et alj|2012l) . How- 
ever, these expectations have been set largely by the assump- 
tion that the stellar mass to halo mass ratio remains a scale- 
free power law below 10 u Mq. As seen in Fig. [JJ the low 
mass power-law behavior is broken below 10 u Mq, corre- 
sponding with an upturn in the stellar mass function below 
10 8 5 M dBaldrv et alj2008h (this result has also been seen by 
Kravtsov, in prep). This underscores the danger of assuming 
that faint dwarfs obey the same physical scaling relations as 
Magellanic Cloud-scale galaxies; moreover, it also is a strong 
argument against fitting the stellar mass - halo mass relation 
with a double power law (see discussion in Appendix ID1. 

Concerning the range of allowed SMHM relations, the ob- 
servational systematics are large enough that our results are 
marginally consistent with an unchanging SMHM relation 
from z = 6 to z = 0. Nonetheless, the feature with strongest 
significance is a gradual decrease in stellar mass in the me- 
dian 10 1 'M q halo from z = to z = 2, followed by an increase 
again for redshifts z > 6; also potentially indicated is an in- 
crease in stellar mass in the median M/, > 10 13 M Q halo from 
z = to z = 2. The best fit SMHM relations at z = 7 and z = 8 
are significantly different than at lower reshifts, with more 



Average Star Formation Histories From Z=0-8 



13 



Time [Gyr] 
4 ' 2 




Time [Gyr] 
4 ' 2 





FIG. 9. — Left panel: Stellar mass growth histories of galaxies of different masses. Lines shown the amount of stellar mass remaining in central galaxies at the 
present day that was in place (in any progenitor galaxy) at a given redshift for our best fit model. Right panel: Amount of stellar mass in the intracluster light 
(ICL) remaining at the present day in place at a given redshift. Note that plots for 10 i5 Mq halos are not shown as they are nearly identical to those for 10 14 Mq 
halos. Shaded regions in both panels shown the one-sigma posterior distribution. 

sistent with all available data. 

5.3. Stellar Mass and Intracluster Light Growth Histories 

Our model constrains the buildup of stars in the intraclus- 
ter light purely from observational galaxy data and measure- 
ments of the halo- halo merger rate in simulations (see also 
IWatson et al.l2012l for an alternate method). In our best-fitting 
model, only 5% of stellar mass in mergers for 10 14 M Q halos 
is allowed to be deposited onto the central galaxy since z = 1, 
and only 10% for 10 13 M© halos. For Milky Way-sized and 
smaller halos, this number rises rapidly to 70-80%. Yet, due 
to the sharply decreasing stellar mass to halo mass ratio for 
lower-mass halos, most of the incoming stellar mass will be in 
(rare) major mergers. Central galaxies are therefore relatively 
uncontaminated by stars from smaller merged satellites. 

In Fig. [9] we show the amount of galaxies' present-day stel- 
lar mass and intracluster light (ICL)/halo stars that was in 
place at a given redshift. The left-hand panel (stellar mass) 
shows that almost all stars in the central galaxy in present-day 
cluster-scale halos were in place at z = 2. However, since that 
time, a large number of their satellites have been disrupted 
into the ICL. Thus, for cluster scale halos, the stellar mass in 
the ICL exceeds the stellar mass in the central galaxy by a fac- 
tor of 4-5, consistent with observations (e.g. jGonzalez et ail 
120051) . We note that our model predicts what may seem to be 
a large ICL fraction for Milky Way-sized galaxies (10 12 M Q ). 
However, the Milky Way is a special case. It has no t had a ma- 
jor merger for ~ 10-11 Gyr (iHammer et al .1120071) : however, 
as noted above, only major mergers can contribute substan- 
tially to the ICL. A major merger 1 1 Gyr ago would have, 
however, contributed less than 3% of the present-day stellar 
mass of the Milky Way into the ICL; allowing for passive 
stellar evolution, this would result in less than 2% of the lumi- 
nosity of the Milky Way coming from the intrahalq light, in 
excellent agreement with observations dPurcell et al.ll2007T) . 

Finally, in Fig. [TO) we show the inferred fraction of stel- 
lar mass growth coming from in situ star formation (as op- 
posed to galaxy-galaxy mergers) as a function of halo mass 
and redshift. At all redshifts greater than 1, the vast majority 
of stellar mass growth is from star formation, and hence stel- 
lar mass functions alone may be used to infer galaxy star for- 
mation rates. However, the stellar mass growth rate of high- 
mass, low-redshift halos is such that it cannot be explained 
entirely by in situ star formation. For these halos, stellar mass 



FIG. 10. — The fraction of stellar mass growth due to in situ star formation 
(as opposed to growth by galaxy-galaxy mergers) as a function of halo mass 
and redshift. 

stellar mass per unit halo mass. However, concerns about the 
r eliab ility of the stellar mass functions at those redshifts (see 
33. U urge caution in interpreting the physical meaning of this 
result. 

A useful perspective on these results can be obtained by 
considering the historical stellar mass to halo mass ratio of 
halos, as shown in Fig. [8] Despite the large systematic un- 
certainties, it is clear that halos go through markedly dif- 
ferent phases of star formation. This evolution is most ap- 
parent for massive halos, as observations have been able to 
probe the properties of the progenitor galaxies all the way to 
z = 8. Specifically, high-redshift progenitors of today's bright- 
est cluster galaxies (M/, ~ 1O 14 M were relatively efficient in 
converting baryons to stars — comparable to the most efficient 
galaxies today. However, between redshifts 2-3, their effi- 
ciencies peaked, and thereafter they began to form stars less 
rapidly than their host halos were accreting dark matter. At 
the present day, such galaxies have an integrated star forma- 
tion efficiency that is two orders of magnitude less than at 
their peak. The picture is less clear for progenitors of lower- 
mass galaxies because current observations cannot probe their 
progenitors as far back. Nonetheless, their apparent behavior 
in Fig.|7]of rising to a peak efficiency and later falling is con- 
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FIG. 1 1 . — Left panel: Ratio of the average SFR to the average baryon accretion rate ( f^dMjdt) for halos as a function of halo mass. Right panel: same, except 
as a function of halo mass at z = (i.e., ratio of SFR to baryon accretion rate for the progenitors of present-day halos). Lines shown the best fit model and shaded 
regions show the one-sigma posterior distribution. 





FIG. 13. — Left panel: Average age of presently-existing stars in halos as a function of stellar mass. Right panel: Time required to form 50% and 90% of 
presently-existing stars as a function of stellar mass. Lines shown the best fit model and shaded regions show the one-sigma posterior distribution. 
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functions alone cannot be used to infer their star formation 
rates (see also Appendix[J]i. Interestingly, in the local universe 
the transition between merger-dominated growth and star for- 
mation growth in the local Universe is roughly the mass of 
the Milky Way — low mass galaxies are dominated by star 
formation, and galaxies in halos more massive than 10 12 are 
dominated by merging. 

5.4. Baryon Conversion Efficiencies 

Assuming that the baryon accretion rate for a halo is equal 
to the universal baryon fraction (fb, which is 0.17 in our cos- 
mology) times the overall halo mass accretion rate, we can 
calculate the instantaneous baryon conversion efficiency im- 
plied by our models. In the left-hand panel of Fig. [IT] we 
show the instantaneous conversion efficiency at fixed halo 
mass. Remarkably, this efficiency is always highest for 
Milky-Way sized halos (10 12 M Q ), corresponding well with 
the peak in the stellar mass to halo mass ratio seen in Fig. 
[7] Moreover, this efficiency (20-40%) is constant to within a 
factor of two over a remarkably large redshift range, suggest- 
ing that the efficiency is only a weak function of accretion 
rate in halos of this mass. For more massive halos, there is 
a decrease in the baryon conversion efficiency following red- 
shift 2-3 by an order of magnitude or more. This suggests 
that the ability of accreted material to cool onto star-forming 
regions becomes dramatically impaired at lower redshifts for 
such halos. However, this decrease is gradual, taking place 
over many Gyr. Rather than an abrupt change in the character 
of infalling gas, this may suggest that the decreasing density 
and accretion rate of gas make it gradually more difficult for 
cold clumps to form, especially in the presence of an active 
galactic nucleus. 

On the right-hand side of Fig. QT| we show the historical 
baryon conversion efficiency for progenitors of z = halos. 
For massive halos at z = 0, their conversion efficiency climbs 
steeply towards higher redshifts as their halo mass falls to 
roughly 10 12 M Q ; once their halo mass drops below that value, 
they then become less and less efficient at very high red- 
shifts. For less-massive halos, such as those that are 10 12 M Q 
at z = 0, their conversion efficiency has been increasing from 
early times to the present day. 

On the left-hand side of Fig.[T2l we show extrapolated his- 
torical stellar mass to halo mass ratios — proportional to in- 
tegrated baryon conversion efficiencies — similar to the bot- 
tom panel of Fig. [7] Notably, stellar mass to halo mass ra- 
tios also peak when halos reach 1O 12 M , indicative of steeply 
falling star formation efficiencies at higher and lower masses. 
It would appear that the maximum integrated stellar mass ef- 
ficiency is around 20-40% at all redshifts. 

5.5. Stellar Ages 

In the right-hand panel of Fig. [12] we show the historical 
stellar mass in progenitors of halos at z = relative to the 
present-day stellar mass. This should not be confused with 
Fig. [9] which shows what amount of the currently-remaining 
stellar mass was in place at a given redshift; these differ 
mainly because of passive stellar evolution (massive stars that 
formed sufficiently long ago will burn out by the present day). 

These data allow us to derive stellar mass-weighted ages 
as a function of stellar mass, shown in the left-hand panel of 
Fig-El] These ages are consistent with the stellar populations 
in massive galaxies (and halos) forming at very early times 
and having little star formation continuing to the present day. 



Less-massive galaxies have younger average stellar ages, con- 
sistent with the ongoing star formation seen in the star forma- 
tion histories for such galaxies in Fig. [6] While we have little 
information on the progenitors of galaxies less massive than 
10 9 M Q , there is evidence that the average ages of the stellar 
populations may increase for such galaxies. Specifically, as 
shown in the right-hand panel of Fig. Q~3] the time required 
to form 50% and 90% of the galaxy's stars increases towards 
lower stellar masses, indicating that the star formation histo- 
ries become more and more fiat (constant), on average. A 
perfectly flat star formation history would result in an average 
age of about half the age of the universe, or roughly 6.7 Gyr; 
this is somewhat higher than that seen for 1O 9 M galaxies. 
On the other hand, for individual dwarf galaxies, the star for- 
mation history is likely to be stochastic, which may result in 
a large scatter around this average value. 

For the most massive galaxies today (~ 1O 12 M ), the time 
required to form 50% of their stars was extremely short — 
only 1-2 Gyr. This in turn translates to a high star for- 
mation rate, on the order of 200-1000 M Q yr , which 
is consistent with observations of ultra-luminous infrared 
galaxies (ULIRGs: iMagnelli et alJ I20TI iDaddi et alJ 120051: 
34J). 



IChapman etail 2004 



5.6. Comparison to Other Results 

We show a comparison of our best-fit results for the stel- 
lar mass to halo mass ratio at z = 0.1 to previously-published 
results in Fig. [14] and we show comparisons at z = 1.0 and 
Z = 3.0 in Fig. Q3] Where possible, conversions to our as- 
sumed cosmology and halo mass definition have been applied. 
The resul ts in this work are alm ost identical to our previous 
results in iBehroozi e t al. (2010), with the exception of a de- 
viation at low halo masses due to the updated stellar mass 
functions used. We also compare to the new constraint of 
iReddick et al.l d2012l) . which uses additional input from the 
correlation function and conditional stellar mass function as 
measured by SDSS. There is a slight discrepancy between the 
two results for small masses due to the different satellite frac- 
tions obtained usin g M nRa k and v„ (still within our system- 
atic errors), but the IReddick et alJ (120121) result is only well 
constrained by additional data in the range 10 12 — — 10 14 . 
Com parisons to other data sets at z < 1.0 dMoster et al 
[GuoeTall I20T0I: IWang & Jingl l20Tot IZhenTetaT 
;t alJ U 




lYang et alJ l2009at lHansen et al.M2009l; iLin & Mohi 
are discussed extensively in IBehroozi et al.l (120101) . so 
we do not repeat that discussion here. The z = 3 comparison is 
notable because it illustrates the large discrepancies that can 
occur when di fferent stell a r mas s functions are used at z > 3. 
For example, lYang et ail ( 120121) perfo rm modeling for two 
separ ate high-z stellar mas s functions dPerez-Gonzalez et alJ 
l2008l and lDrorv et alj|2005l) . Part of the discrepancy between 
the two sets of results may be due to a somewhat restrictive 
redshift fit, and part of it may be due to systematic biases in 
stellar masses at high redshifts. 

In Fig. 021 we show a comparison of galaxy st e llar m ass his- 
tories between our results and thos e of | Leitnerl d2012l) . A di- 
rect comparison is difficult because ILeitnerl (120121) only con- 
siders star-forming galaxies. However, according to Eq. [8] 
the active fraction for 10 9 5 M Q galaxies is 90%, and the ac- 
tive fraction for 1O 1O M galaxies is 65% — although most of 
the passive 1O 1O M galaxies have only recently become pas- 
sive. Comparison of the stellar mass histories for those galax- 
ies is thus somewhat more feasible. Our stellar mass histories 
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FIG . 14. — Comparison of our b est-fit model at z = 0.1 to previously publ i shed results. Results compared include those from our pr evious work IBehroozi et al.l 
1201(1). from abundance matching Moster et al.ll2l)12r . IReddick et al.]|2012tlMoster et alj|2(JT0l:lGu o et al. 2010; Wang & Jing 20J3), from HOD/CLF modeling 
IZheng et al. 2007a; Yang et al. 20 1 2 ), and from clust er catalogs (Yang et al. 2009a; Hansen et al. 2009; Lin & Mohr 2004j). Grey shaded regions correspond to 
the 68% confidence contours of Behroozi et al. (2010). The one-sigma posterior distribution for our model is shown by the red error bars. 
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FIG . 15. — Comparison of our best- fit model at z = 1.0 and z = 3.0 to prev i ously published results. Results compared in clude those from our previous 
work (B ehroozi et al. 20i rj), from abundance matching (Moster et al. 20ia i20Trj| ; [Conrov & Wechsleil20 09; Wang & Jing|2010l) . and from HOD/CLF modeling 
I Zheng et al. 2007a; Yang et al. 2012; Wake et al. 201 1 ). Yang et al. (2012) reports best fits for two separate stellar mass functions, and we include both at z, = 3.0. 
Grey shaded regions correspond to the 68% confidence contours ofjBehroozi et al. 1 2 0101) . 
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FIG. 16. — Com parison of our best-fit model's stellar mass histories to those 
of lLeitneil 120121) (lines; shaded reg ions show the one-sigma posterior distri- 
bution of our model). Note that theJLeitner (2012) results only consider star 
forming galaxies; ours show the average histories for all galaxies. 




< 0.86 



FIG. 17. — Allowed range for the completeness (i.e., probability of galaxy 
detection) of the surveys we use as a function of redshift; the black line shows 
the median value, and the grey bands show 68% confidence contours. This 
average applies only to galaxies above authors' stated completeness limits. 
No strong evidence is seen for missing large numbers of high-redshift galax- 
ies due to burstiness or dustiness. 



agree remarkably well with those in iLeitnerl (120 121) except at 
recent times for 10 10 M Q galaxies (where our model includes 
more passive galaxies) and at early times for 10 9 5 M Q galax- 
ies, where observations do not constrain properties of their 
progenitors. 

5.7. Systematic Uncertainties 



Several aspects of the allowed parameter range are notable. 
Most importantly, there is no necessity for a large number 
of galaxies to be missed (due to dustiness or burstiness) at 
high redshift; indeed, the allowed incompleteness is in gen- 
eral less than 30%, as shown in the top panel of Fig. [17] This 
would imply that most high-redshift surveys are not missing a 
large fraction (>50%) of galaxies above their nominal detec- 
tion thresholds. 

The allowed parameter range is consistent with having both 
[i = and k = at z = 0; i.e., no systematic offsets in stellar 
mass are necessary at z = 0. This is partially by design: the 
parametrization of the SMHM was chosen so that it could fit 
the z = constraints without need for fi or k; moreover, the 
cosmic and specific star formation rates are derivative con- 
straints, meaning that they only constrain how the SMHM 
evolves with time. While the main fits favor some redshift 
evolution in the systematic parameters, it is still possible to 
reasonably fit all observations while disallowing all system- 
atic corrections (i.e., assuming that the stellar masses and star 
formation rates used in this paper are on average correct and 
that all the surveys were 100% complete). Details of this 
fit are presented in Appendix |J] We also present discussions 
therein of how the allowed parameter space varies when cer- 
tain data sets are not used to constrain the fit; specifically, 
excluding z > 6 data; excluding cosmic star formation con- 
straints; and excluding all star formation constraints. 

5.8. Fits to Star Formation Histories 
The most commonly used star formation history fitting for- 
mula is a declining exponential with time (SFH(t) oc e~'/ T ), 
used by almost all the observational papers cited in this work 
(see references in Tables [3] and |6j. Various modifications 
to this form have b een suggested (e.g., te~'l T or e + 'l T in 
iMaraston et al.ll2010l) . and there has been increasing recent 
evidence that increasing star formation histories are more ap- 
propriate for massive galaxies at high redshifts. 

Our best constraints in this study are on average star forma- 
tion histories for halos at a given mass. Thus, we can place 
the strongest priors on the expected form of stacked galaxy 
star formation histories. However, we can also place some 
weaker constraints on what individual galaxy star formation 
histories should look like. 

Individual galaxies' histories can differ from the average 
star formation histories for three reasons. The first is stochas- 
ticity in the star formation rate on short timescales relative 
to the dynamical time of the halo, due to feedback effects as 
well as shot noise in the star formation rate. This variation 
cannot be modeled with our current approach, so we caution 
that the constraints on individual galaxy histories that we de- 
rive should be considered as smoothed over a period of at least 
500 Myr. If enough color bands are available to avoid degen- 
eracies with other fitting parameters, observers may consider 
adding the possibility of recent starbursts to their models to 
account for the fact that they can significantly alter observed 
galaxy colors. 

The next reason for deviations is sustained deficits in the 
star formatio n rate for satellite galaxies compared to field 
galaxies ( e.g. IWetzel et al.ll2012l) . As shown in Fig. [24] (see 
also e.g. iRravtsov et al.l 12004b . the satellite fraction of ha- 
los falls both with increasing mass and increasing redshift. 
This is not a practical concern for z > 3, where the satellite 
fraction is less than 20% for all currently-observable galax- 
ies (M„ > 10 9,5 M Q ). At lower redshifts, the observation that 
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FIG. 18. — Featureless, rising power law star formation histories are appropriate for z > 3; however, for z. < 3, galaxies have peaks in their star formation rates 
that depend on halo mass. These figures show constraints on individual star formation histories of galaxies in halos of mass 10 11 - 1O 14 M as a function of 
time from the beginning of the universe to z = 0.1. These halo masses correspond to stellar masses of ~ 10 9 , 3 x 1O 1O M , 10 n M o , and2x 10 n M Q aU = 0, 
respectively. 
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FIG. 19. — Left panel: Median mass accretion histories for halos in narrow mass bins (0.25 dex). The shaded regions contain 68% of the spread in mass 
accreti on histories for different halos. Right panel: Constraints on individual star formation histories for galaxies in 10 12 Mq halos, compared with fits from Eq. 
1221241 

satellites cease to fo rm stars relatively quickly after accretion 
(IWetzel et alJl2012h means that declining star formation his- 
tories may continue to be reasonable fits for satellite galaxies, 
even though most other galaxies of the same mass may be 
better fit by rising star formation histories. 

Finally, we note that differences can arise in individual 
galaxy star formation histories because of the scatter in stel- 
lar masses at a given halo mass, as well as the scatter in mass 
accretion histories for halos (see Fig. [19}. We can attempt to 
model these effects by sampling random mass accretion his- 
tories from the Bolshoi simulation and then sampling random 
stellar mass growth histories as allowed by the scatter in stel- 
lar mass at fixed halo mass (i.e., choosing random stellar mass 
offsets at z = 8, z = 1 .0, and z = and using spline interpolation 
at intermediate times). Our results, expressed as a function of 
time since the beginning of the universe, are shown in Fig. [18] 

These results have dramatically larger error bars than those 
in Fig.[6]on account of the spread in mass accretion histories 
for halos. Nonetheless, some basic conclusions can be drawn. 
Star formation histories for z > 3 galaxies increase with time. 
Although there are different ways to parametrize these histo- 
ries, a straightforward one is a direct power law form: 



of iPapovich et al.l (1201 11) ; however, we find steeper slopes 
(B ~ 3 - 4) than theirs (B ~ 1 .7) because they ignore all ef- 
fects of mergers that occur from z = 8 to z = 3. 

At lower redshifts, there is a mass-dependent turnover af- 
ter which the star formation history begins to decline. This 
happens at z ~ 2-3 for lO n M galaxies, making a declining 
exponential at z = a reasonable fit. However, it happens as 
late as z = 0.5 for 10 9,5 M Q galaxies — meaning that a declining 
exponential is never a good fit for these galaxies unless they 
are satellites. The best fit in general for all the constraints on 
individual histories presented here is a double power law: 



SFH{t)=A 



(23) 



typical values of Z? and C range from 1 to 5 . It may often be the 
case that the available data is insufficient to constrain all three 
shape parameters. In this case, a hybrid of the exponential 
decline and power law rise still provides a reasonable fit: 



SFH(t) = At B exp(-t/T). 

6. DISCUSSION 



(24) 



SFH(t) = At (forz>3). 



(22) 



This is also equally capable of fitting the average star for- 
mation rates in Fig. [6] at z > 3. Our result is similar to that 



The existence of a "cold mode" of gas accretion that al- 
lows efficient star formation at high redshifts and shuts off for 
massive galaxies past z ~ 2 has been pre dicted in hydrody- 
namical simulations for the past decade dBirnboim & Dekel 
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[2003b . In this work, we observe a transition after z ~ 3, where 
galaxies in massive halos become progressively less efficient 
at forming stars, even after correcting for slowing accretion 
rates since that time (Fig. Q~T] right panel). At fixed halo mass 
(Fig. [TT] left panel), the transition is less severe, but there 
has been a clear reduction in the conversion efficiency since 
z ~ 3 to the present day. This suggests a picture in which 
dense infalling gas can cool to the galaxy disk efficiently at 
high redshifts; however, at lower redshifts, the infalling gas 
becomes sparse enough that it is susceptible either to shock 
heating or feedback from an active galactic nucleus. 

Another important feature of our results is that merging 
galaxies' stars almost never reach the central galaxy in mas- 
sive halos (< 10% si nce z = 1 for halos larger than 10 13 Mf^) . 
This claim is not new dPurcell et al.l2007tlConrov et al.120071) . 
but it is reassuring that we duplicate this conclusion with 
hig h confidence. This is in stark contrast to what is asserted 
in dMoster et al.M2.Q12b (i.e., 80% of merging stars reach the 
central galaxy); however, their assertion is based on sim- 
ulations of Milky-Way sized halos, for which we find that 
indeed ~ 70% of merging galaxies' stars reach the central 
galaxy. Thus, the fraction of merging stars that reach the cen- 
tral galaxy is likely a strong function of halo mass. 

The amount of stellar mass that was contained in disrupted 
galaxies can easily reach five times the amount of stellar mass 
in the central galax y for high-mass halos (Fig. |9), consistent 
with observations in lGonzalez et al.l (120051) . Regardless of un- 
certainties in determining stellar masses, there is no possible 
way to funnel any significant fraction of this mass into the 
central galaxy without badly violating observed limits on the 
number density of massive galaxies. Rather, stars from merg- 
ing galaxies must end up in a more extended distribution to 
avoid confusion with the central galaxy, a topic that will be 
explored in depth in a future paper (Behroozi et al., in prep). 
From the opposite perspective, galaxy stellar populations con- 
tain mostly stars from the most-massive progenitor, with some 
contamination from major mergers. Surveys of stars ithin a 
galaxy's optical radius (such as for the Milky Way) thus cap- 
ture a potentially limited range of different star-forming envi- 
ronments. A wider range of environments can be sampled by 
considering the ICL and halo stars in general. 

The most efficient halo mass for converting baryons into 
stars appears to be close to 10 12 M Q , regardless of redshift. 
The progenitors of today's massive halos (> 10 I4 M Q ) passed 
through this mass quickly, and as a result formed most of their 
stars in a short time: a mere 1-2 Gyr (Fig. fl3T) . Given that 
baryon to star conversion becomes so inefficient above and 
below this halo mass, much of the shape of the star forma- 
tion histories seen in this paper (Fig.(6jl could be explained by 
considering the duration of time that a given halo's mass ac- 
cretion history spends close to this most efficient halo mass. 
That is to say, the peakiness of the star formation history de- 
pends on how quickly the halo mass accretion history passes 
through halos of mass ~ 10 12 M Q . As lower -mass halos have 
slower relative accretion rates, they have flatter star formation 
histories; e.g., Milky Way-sized halos have lingered in this 
mass range from z ~ 1, and so have the current highest stellar 
mass to halo mass ratio, as well as star formation histories that 
track their halo mass accretion rates well at late times. 

Low-mass halos (10 n M©) show an unusual behavior in 
terms of their baryon conversion efficiency. In the left panel 
of Fig. QT] it would appear that the instantaneous baryon con- 
version efficiency of 1O U M0 halos is high at high redshifts 



(z > 5), low at intermediate redshifts (2 < z < 5), and then 
high again at low redshifts. While this behavior matches what 
is seen in the integrated baryon conversion efficiency for such 
halos (Fig. [12] left panel), the corresponding galaxies are of- 
ten at the faint edge of the completeness limit for stellar mass 
surveys. For that reason, observational biases may be the most 
likely source of this behavior. 

Our results have importance for future observational stud- 
ies. The current levels of systematic errors in the inference of 
stellar masses and star formation rates are the dominant sys- 
tematic uncertainty in this work. As discussed in §0 stellar 
mass functions at z > 3 can disagree by up to 0.5 dex, well be- 
yond any errors expected from sample variance. Star forma- 
tion rates are equally concerning, both at high redshifts and 
at low redshifts, especially in terms of specific star formation 
rates. As a result, systematic errors are the single most impor- 
tant aspect preventing better understanding of t he f ormation 
histories of massive galaxies. As discussed in §15.81 we have 
provided new fitting formulas for galaxy star formation his- 
tories, which should help reduce one of the largest sources 
of systematic error especially for low-mass and high-redshift 
galaxies. New studies of stellar mass functions and star for- 
mation ratesjromjle^gnnj^^ CAN- 
DELS dKoekemoer et al.ll201 U [Grogin et al.ll201 ll) . analyzed 
with self-consistent star formation histories and assumptions 
about systematic errors, will provide a large step forward in 
constraining the model presented here. 

7. CONCLUSIONS 

We have presented a comprehensive analysis of average star 
formation rates and histories in galaxies from z = to z = 8, 
and their connection to the underlying growth and merging 
of dark matter halos, along with a treatment of the inherent 
uncertainties. Our approach provides a self-consisent picture 
of the growth of galaxies over this epoch in the cosmological 
context of ACDM structure formation. The model is able to 
match the evolution of cosmic star formation rates, specific 
star formation rates, and the stellar mass function of galax- 
ies over the last 13.2 billion years of cosmic time. Our main 
findings are as follows: 

1 . Halos of mass <~ 10 12 M Q appear to be the most efficient 
at forming stars at every epoch, with baryon conversion 
efficiencies of 20-40% over nearly the entire redshift 
range of this study. Halos at higher and lower masses 
are less efficient by orders of magnitude, especially at 
low redshifts. 

2. The baryon conversion efficiency of more massive ha- 
los is still reasonably high (10-20%) until z = 2-3, after 
which time it takes a steep downturn. This is consistent 
with expectations of reduced cold mode accretion onto 
the corresponding galaxies. 

3. We have characterized the fraction of stars in galax- 
ies due to merging vs star formation. At z > 1, 
galaxy buildup in all but the most massive galaxies was 
strongly dominated by star formation. At the present 
day, the transition between merger-dominated buildup 
(at high mass) and star formation-dominated buildup (at 
low mass) occurs at roughly the mass of the Milky Way. 

4. We confirm previous results in support of most merging 
galaxies' stars being disrupted in the ICL for massive 
halos. We find, however, that the fraction of merging 
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stars that reach the central galaxy may increase strongly 
with decreasing halo mass, to about 70% for Milky 
Way-sized halos. 

5. For massive galaxies (M* > lO n M ), we predict that 
at least as much stellar mass is present in the ICL as in 
the galaxy itself out to intermediate redshifts, possibly 
as high as z ~ 1 for the largest clusters (although this 
depends on the exact cut between ICL and galaxy stellar 
mass). 

6. Star formation histories peak at z ~ 3 for massive galax- 
ies and are increasing nearly to the present day for small 
galaxies. We have presented new fitting formulas for 
galaxy star formation histories: a featureless power law 
for high redshifts (Eq. [22b . and more general formulas 
valid for all redshifts (Eq.|23]and|24j. 

7. Upturns in the stellar mass function at low stellar 
masses mean that dwarf galaxies (M* < 10 8,5 M Q ) have 
a higher stellar mass to halo mass ratio than that 
predicted by extrapolating results from more massive 
galaxies. The upturn is most robustly observed at low 
redshifts, but there are also indications at higher red- 
shifts for a steepening slope to the SMF at low stellar 
masses. The stellar mass to halo mass relation thus can- 
not be fit with a double power law alone. 

8. Systematic errors remain at high levels (~ 0.2-0.3 dex 
in stellar masses, 0.15-0.3 dex in SFRs), and disagree- 
ment between galaxy abundances for high-redshift sur- 
veys (z > 3) is substantial. 

9. Galaxies with stellar masses below 10 10 M Q have the 
least well-constrained properties in this study. This 



could be improved by future surveys and by incorporat- 
ing information on the star formation histories of local 
galaxies. 
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APPENDIX 

IMPACT OF ASSUMING AVERAGE STAR FORMATION HISTORIES ON DERIVED STAR FORMATION RATES 

As mentioned in the introduction to ^2] there is a potential error in our method because we propagate average star formation 
histories for halos at a given mass and redshift, as opposed to individual star formation histories. We discuss effects from both a 
distribution in star formation rates for galaxies at a given halo mass and from a distribution of halo mass accretion histories. 

Bimodality has been observed in the s tar formation rates of galaxies at fixed stellar mass, partially but not exclusively driven 
by quenched satellite gala xies (see, e.g.. lWetzel et alJl2012l) . This may appear to be a problem for our method because the rate 
of specific mass loss (see 32.3b will be different depending on whether a given galaxy is active or passive. Remarkably, however, 
our method is robust to any distribution of star formation rates in galaxies of the same stellar mass. Because our method only 
depends on the total stellar mass buildup in halos of a given mass Mi, at any given time, the fact that some galaxies get a larger 
share of that buildup than others is not a concern. If certain galaxies sustain a higher growth rate for a long time, then the scatter 
in stellar mass at fixed halo mass will certainly increase; however, the possibility of scatter changing with time is an effect that 
we include explicitly in our model. Bimodality in star formation rates is only a problem if it results in a distribution of stellar 
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masses at fixed halo mass that cannot be reasonably modeled by a log-normal distribution (i.e., the d istribu tion we assume in this 
work). So far, no evidence has emerged for any other distribution to be preferred (see comments in $2.2. II where £ is defined). 

There is more possibility for error because we do not account for correlations between the star formation history (SFH) and 
the mass accretion history of a halo. For example, a rapidly accreting halo might have a more recent star formation history than 
oth er halos at the same mass. In this case, when it transitions from one mass bin to another, we may use the incorrect SFH in Eq. 
IB2I This effect is balanced by halos that accrete more slowly than other halos of the same mass, which cancels the effect to first 
order. 

It is worth working through a more quantitative example to verify this assertion. Consider a population of galaxies all with 
the same stellar mass M* = 10 Mq today. Assume that each galaxy has had a constant SSFR for its entire lifetime (i.e., an 
exponential growth history), and moreover that the distribution of the SSFRs is log-normal with width 0.3 dex (similar to that 
seen in lSalim et al.l |2007l) and median 10~ 9 yr" 1 . The average (as opposed to median) star formation rate for this population is 
12.7 M Q yr" 1 ; given the star formation histories just specified, this translates into an average rate of change of stellar mass of 
9.3 Mq yr" 1 over the next 100 Myr. If, on the other hand, we had simply assumed that all of the galaxies had star formation 
histories corresponding to a constant SSFR of 10 yr" 1 , we would interpret a buildup in stellar mass of 9.3 Mq yr" 1 as implying 
an instantaneous star formation rate of 13.4 Mq yr" 1 , an error of 6%. Thus, we conclude that the errors on the average quantities 
considered here introduced by lumping together star formation histories for halos with different mass accretion histories are 
strongly subdominant to existing observational uncertainties. 

MERGING GALAXIES AND MASS ACCRETION 

Because of passive stellar evolution, calculating the amount of stellar mass remaining in a galaxy (required for determining 
the star formation rate in Eq. [D, one must know the entire star formation history. For a large simulation, keeping track of the 
individual stellar mass histories (covering hundreds of timesteps) for millions of galaxies repeated over several million realizations 
of the stellar mass - halo mass relation is inefficient. Instead, since we are interested in only the average star formation rate as a 
function of mass, we take an elegant (although approximate) shortcut. 

Instead of keeping track of every halo/galaxy individually, we instead keep track of average star formation histories at each 
timestep in bins of halo mass (which is a rank-3 tensor, SFH'/°""; t is the time now, tf orm is the timestep at which the stellar 
mass was formed, and m is the halo mass). For central halos, we use the virial mass (as discussed in §|4); for satellite halos, we 
use the mass of the largest progenitor (the peak mass) in the satellite's mass accretion history. At every timestep, we calculate 
the number of halos that transition between one mass bin and another (T t l^), i.e., the number of halos that had a mass m p at 
timestep t and a mass m„ at timestep t+ 1 (see Appendices|G]and|H|. Similarly, we calculate the merger rates for subhalos at each 
timestep as a function of the subhalo peak mass and the target host halo mass (M^L), as discussed in AppendixU Keeping track 
of star formation histories then becomes an exercise in tensor multiplication, which is fast regardless of the number of halos in 
the simulation. 

To convert Eq. [T]into a full tensor equation, we need a small amount of additional notation. The stellar mass at a particular 
timestep t and halo mass m becomes SM U „, the star formation rate at a particular timestep and halo mass becomes SFR, m , and 
the number counts at a given timestep for halos in a given mass bin becomes N ljn . In addition, we denote the difference in stellar 
mass from the previous timestep by ASM ljn . Thus, Eq.Q]formally becomes: 

ASM,. m = ^SM,. m - J2 SFH'/;r ( 1 - Ci oss (/, t form - At , t form )) 

SFR f S FR(m,t)ASM,, m 

Af(l-C loss (f,f-Af,f)) 

Closs(^«ow it a i tb) — \tb ta\ I /loss (jnow X)dx^ (Bl) 

where Q s S (t now ,t a ,ti,) is the average stellar mass loss fraction (a.k.a., "recycling" fraction) at the present time for stars that formed 
at a constant rate from t a to (Note that the indices t and m are not meant to be summed over). 

To calculate st ar fo rmation histories, we need to know both internal star formation and the amount of stars accreted in mergers. 
As discussed in 32.41 purely definitional changes in what constitutes the boundary between a galaxy and the ICL can effectively 
shift stellar mass from the ICL into the galaxy. Not knowing the relative fraction of each, we assume that merging stellar mass 
from satellites takes priority over stellar mass from the ICL. For convenience, we define SM' t m and ICL 1 , m be the total stellar 

mass in merging satellites and the total stellar mass in the main progenitors' ICL, respectively; we also define ICLH'/""" to be the 

star formation history of the ICL (analogously to SFH'/f,""). The new star formation histories will then be the sum of the old star 
formation histories from all progenitors, with the addition of new stars formed at the latest timestep and material from mergers 
and the ICL: 

sFH;ij lm =Y,[sFH!x'{T,7+cr;:' se M m) 

m p 

emerge ■ f, ( 1 ~ fsFIt(m, t )) ASM t+Al . m \ 

C <>>» =mm ( v 1 ' SMl m J (B3) 
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rdef f n ( 1 ~ fsFRjm, t))ASM l+Aljn - SM' tM \ 

C,J = max [0, ^- j (B4) 

SM' t m = SFH'/^MZ^ ~ QoJf + At , t form - At , f /o ,,„)) (M +A ,, m )~ 1 (B5) 

m p 

ICL' l m = ^ ICLH'/Z* T™£ ( 1 - Cioss(* + Af , t form - At , t form )){N, + A,,mT l • (B6) 

(Again, the indices t and m are not meant to be summed over; 8/""" is the standard identity tensor: 8/° m = 1 if t = tf orm and 
otherwise). Not e that, in this present analysis, we do not separate the star formation histories of satellites and centrals. As 
discussed in 35.81 satellites at high redshifts (z > 1) tend to have short infall timescales, so they have roughly similar star 
formation histories as centrals. At later times, this approximation is less valid; however, we can nonetheless correctly recover the 
average star formation history across all halos (centrals and satellites) at a given halo mass, as discussed in AppendixlAl 
Calculating the ICL formation history as well as the total amount of stars in the ICL is similar: 

ICLH'/ + °Z t , m = Yl (sFHl%(l-C^lift+ICIJti%T%(l-C%)) {N t+At ^f x (B7) 

m p 

ICL t m = 2^ ICLH'j^" ( 1 - Cioss(* ) tform ~ At , tform)) ■ (B8) 

tform 

Note that the merg ing s atellite ICL i s not deposited into the central galaxy's ICL, which is the same approach taken by 
IConrov et al.l (120071) and iPurcell et al.l (120071) . Because the ICL of a merging satellite halo will get stripped long before the 
stars in its associated galaxy, the stripped satellite's ICL will end up in a much more extended, faint distribution than the host's 
ICL. 

This formal specification is equivalent to a discrete differential equation for the stellar mass history. It is thus important to 
properly specify the boundary conditions; that is, the stellar mass history at t = 0. In the absence of constraints on stellar mass 
functions above z ~ 9, we extrapolate the fit to the stellar mass - halo mass relation to z = 15 and assume a uniform star formation 
rate prior to that redshift. While incorrect in detail at z = 15, this choice has little effect on the star formation rates or histories of 
galaxies at later (and more visible) times. 

CALCULATING OBSERVABLES 

Given the large number of equations in our framework, it is appropriate to describe exactly how a point in the parameter space 
of the cosmic stellar mass history (Eqs. |3]|5]l and uncertainties (Eqs. l9lfT8l gets converted into quantities comparable to published 
observations. The three observables in our case are stellar mass functions, specific star formation rates, and the cosmic star 
formation rate density. We deal with each in turn. 

The Observed Stellar Mass Function 

Equation [3] gives the median stellar mass as a function of halo mass for the entire redshift range we consider. If no scatter 
between stellar mass and halo mass were present, the conversion from the halo mass function to the fiducial stellar mass function 
(4>fid) would be straightforward application of the chain rule: 

, , N dN dN ( d\og m M.,\~ X 

<^d(M*)=- — = - —— - 810 . (CI) 

dlog 10 M* d\og 10 M h (M*) \dlog 10 M h J 

In other wo rds, the stellar mass fu nction is equal to the halo mass function divided by the logarithmic derivative of Eq. [3] As 
discussed in lBehroozi et al.l (120101) . this may be convolved with the scatter in observed stellar mass at fixed halo mass (P(Am„ )) 
to give the true stellar mass function: 



<?W(M*)= / fid (M„ + A M JP(A M ,MA M ,. (C2) 

J-M, 

In our model, both the scatter in true stellar mass at fixed halo mass (Eq. [5]) and the scatter in observed stellar mass at fixed 
true stellar mass (Eq.fTTI) are log-normal distributions; hence P(Am,) is a log-normal distribution with width equal to the sum in 
quadrature of the widths of the two sources of scatter ( y / a 2 + £ l 2 ). 

To convert from the true stellar mass function to the observed stellar mass function, we need to correct not only for observational 
completeness (Eq. [T3l). but also for systematic biases in the measured stellar masses (Eqs. l9land[T0T). This latter correction is 
made somewhat trickier by the fact that we only know the passive fraction as a function of the measured stellar mass. The true 
passive fraction, /passive.trueC^*) is related to the observed passive fraction by: 

f {A/T \— 0true(^*,meas 10 ^)./passive,true(^*,meas 10 ^) 

/passi *' = 0true(M t , me a S lO-^)/ p a ss ive,true(M,, me a S lO^) + tme (M,, me a S lO-^--)(l " /pa S .sive,true(M„, meas 10"^)) ' 

In the case that k = 0, then this equation reduces to / P assive,trae(A^*,meas lO - ' 1 ) = / P assive,obs(A^*,meas)- Otherwise, this equation becomes 
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a recurrence relation: 

[(1 -/passive.obsCM*,^))" 1 - l] . (C4) 



/passive, true(-^*,meas 10 ^) _ 0trueC^*,measlO ^ ) 



1 /passive, true(^*,meas :ik'.i- 

10"^) 

Solving this equation is possible, but using the result has disadvantages. Specifically, calculating /passive.mie couples the values of 
ill-constrained regions of the stellar mass function to well-constrained regions through this recurrence relation. The values for 
a — the faint-end slope of the SMHM relation-and n will then be constrained by the MCMC algorithm attempting to improve fits 
to the well-constrained regions of the stellar mass function, and will lose their physical significance. 

We therefore compromise by decoupling the true passive fraction from the value of k. We take /p a ssive,tiue(^*,measlO~ AI ) = 
/passive,obs(M H , iin eas), the exact solution for k = 0. This results in a change to the definition of k so that it no longer corresponds 
to an exact offset for active galaxy stellar masses; however, the correlation between k and the average offset implied for active 
galaxies should still be high. With this approach, it is straightforward to calculate the measured stellar mass function: 

^measC^'K.meas) = /passive, true(^*,meas 10 ^)0true(^ / ^*.meas 10 ^) ~t~ ^trueC^+.meas 10 ^ )(1 /passive. true(^*,measl0 ^ )). (C5) 

Or, in terms of the observed passive fraction, 

^measC^+.meas) — /passive, obsC^*,meas) ( ArueC^*,meas 10 ) + <^trueC^*.meas 10 ^ )(1 /passive, obsC^*,meas 10 )). (C6) 

Then, to correct for observational completeness, we multiply by the appropriate factor: 

0obs(M*,meas) = c(z)4>meas(M*, me - ds ) . (C7) 

Finally, as surveys usually return the stellar mass function over a range of redshifts (e.g., zi < z < Zi), we must perform one 
last integration to determine the survey mass function: 

*,meas )dV c (z) 

0survey(^*,meas) — 7T~, 7 777 T ; (C8) 

V c (Z2)-V c (zi) 

where V c (z) is the comoving volume out to redshift z- 

Calculating the Observed Specific Star Formation Rates and Cosmic Star Formation Rate Density 

Eq. IB II gives the average star formation rate for halos as a function of halo mass and redshift. Eq. [3] gives the median stellar 
mass for halos as a function of halo mass and redshift. Ignoring scatter, the fiducial SSFR is simply: 

SSFR M( M.^ SFR Z h(M - )) . <C9) 

Scatter in the observed stellar mass at fixed halo mass results in contributions from lower- and higher-mass halos, weighted by 
number density. If there is no correlation between stellar mass and star formation rate beyond the correlation between halo mass 
and star formation rate, then the specific star formation rate would be as follows: 

/oo 
(MM, + A M , )SFR{M,,(M* + A„, ))P(A„. )dA M , . (C10) 

However, with a nonzero covariance (p) between SFR and stellar mass at fixed halo mass this equation is amended to become: 

/oo 
(MM. + Am, )(SFR(M h (M* + A„, ))lQ pA '^ )P(A M , )dA M , 
M, 

A ^ =1Og H M.exp(0.5(einl0) 2 ) )- (CU) 

The formula for A' u is necessary because M/,(M» + Am») gives the halo mass at which the median stellar mass is M* + Am„; 
however, the covariance between SFR and SM is calculated relative to the average stellar mass at fixed halo mass. The exponential 
in the denominator (exp(0.5(£ln 10) 2 )) is exactly the ratio of the average value to the median value of a log-normal distribution 
with standard deviation £. 

There are restrictions on the plausi ble range of covari ance coefficient p. At low redshifts, a multiple regression analysis on 
group catalogue data from the SDSS ( Wet zel et al.ll2012h suggests that, at fixed halo mass, the average star formation rate rises 
by 0.23 dex per 1 dex increase in stellar mass. At high redshifts (e.g., z > 6), the covariance will be almost perfect (1 dex rise in 
SFR per 1 dex rise in stellar mass) because of the short times that galaxies have had to form stars. Between these two times, it 
is reasonable to assume that the covariance increases with increasing redshift, so we adopt the following formula to interpolate 
between z ~ and z > 6: 

p(a)= l + (4po.5-3.23)a + (2.46-4p .5)fl 2 - (C12) 

This formula for p(a) gives p(0) = 1, p(l) = 0.23, and p(0.5) = po.5, where po.5 is a free parameter in our analysis. 

By comparison, the total (cosmic) star formation rate density is straightforward to calculate; it is simply the sum of all star 
formation in all halos: 

CSFR lme = [ SFR(M h ) f N d\ogM h . (C13) 
J.^ dlogM h 
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FIG. 20. — Left panel: The stellar mass-halo mass relation as obtained by deconvolution of observed data iMoustakas et al. 2012; Baldrv et al. 2008 ) with an 
assumed scatter in observed stellar mass at fixed halo mass of 0.2 dex; the dashed and dashed-dotted lines show a double power law fit and the five-parameter fit 
used in this work. The double power law fit differs by up to 0.1 dex from the observationally-derived relation. The effects on the stellar mass function are shown 
in the right panel; namely, the fits shown in the left panel have been convolved with 0.2 dex scatter as well as with the halo mass function so as to compare to the 
original stellar mass functions. Again, differences of 0. 1 dex are seen in the double power law case, which is well outside of the reported error bars. 

To calculate the measured quantities (as opposed to true quantities), we adopt the same conversion between measured stellar 
mass and true stellar mass as for the stellar mass function. However, we do not impose additional corrections for active galaxies; 
by definition, because recent star formation is more evident in such galaxies, it is somewhat easier to recover the star formation 
rate, as opposed to the overall stellar mass. Therefore, we only apply a uniform systematic correction and do so in the same way 
as for stellar masses: 

SFR m UM h ) = SFR^W. (C14) 



To calculate SSFR meds and CSFR meds , we would therefore replace SFR tme with SFR meds in Eqs. lCl ll and lC13l There is one more 
complicating factor; namely the burstiness / dustiness of high-redshift galaxies. To account for this, we add additional factors to 
calculate the observed SSFR and CSFR: 



SSFR obs (M*) = 



1 



-SSFR^M.) 



l-b(l-c(z)) 
CSFR obs = [1 - (1 -b)(l - c(z))]CSFR D 



(C15) 



THE USE OF DOUBLE POWER LAWS TO FIT THE STELLAR MASS - HALO MASS RELATION 

A double power law has been used to fit the stellar mass (or luminosity) to halo mass relation at least since IYangetal.l(l200l . 
While using it is simple and convenient, the quantity and quality of galaxy observations have improved dramatically since its first 
use. To test its continued usefulness, we have used abu ndance matching to extract the implied stel lar mass - halo mass (SMHM) 
relation from our low-redshift stellar mass functions (Moustak as et alj|2012t iBaldrv et al.1 12008). To account for the presence 
of scatter in observed stellar mass at fixed halo mass, we use a method based on Richardson-Lucy deconvolution to extract the 
median stellar mass as a function of halo mass. We assume a combined scatter of 0.2 dex. 

The deconvolved SMHM relation is shown in Fig. [20] The figure demonstrates that a double power-law can approximate the 
relation, but can be discrepant by up to 0. 1 dex. Adding one more parameter improves the fit substantially, to a mean offset 
of less than 0.025 dex in our case. The discrepancies in the SMHM relation translate directly to discrepancies in the assumed 
stellar mass function. This may be seen in the right-hand panel of Fig. [20j the fits for the SMHM relations have been used to 
populate the halo mass function and reconvolved with 0.2 dex scatter to determine the resulting stellar mass functions. While the 
five-parameter fit that we use is largely within the error bars of the observational data, the double power law fit gives results that 
are again discrepant by up to 0. 1 dex or 4<r from the observed stellar mass function. 

It is tempting to ignore this issue: after all, the size of the systematic uncertainties in stellar mass is significantly larger (^0.25 
dex). However, if it is indeed the case that the deviations from the power law are driven by systematic uncertainties, then part of 
the value of the physical interpretation of the double power law is lost: the high-mass and low-mass slopes are now contaminated 
by using the incorrect parametrization. In our five-parameter fit, the high-mass dependence and low-mass slope are somewhat 
more insulated from unparametrized systematic errors in the transition region, leading to a more fair estimation of the high- and 
low-mass behavior of the SMHM relation. 

Naturally, no fit will perfectly reproduce the original function. Special caution is necessary when using a MCMC method where 
a parametrization to tightly-constrained data is accompanied by additional systematic or nuisance parameters. This is because 
the MCMC method may try to recruit the additional parameters to improve the fit, rather than for their intended purpose (i.e., 
to cover the range of allowed uncertainties). As shown in Fig. [20] our chosen parametrization can reproduce the local SMF to 
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M h [M e ] MJM Q ] 

FIG. 2 1 . — Left panel: the stellar mass to halo mass ratio inferred via direct abundance matching of stellar mass functions (Table[5J to halo mass functions from 
z = to z = 8, with zero scatter. Right panel: same, with deconvolution of 0.16 dex intrinsic scatter (in stellar mass at fixed halo mass) and (0.07 + 0.04z) dex of 
measurement scatter (in measured stellar mass at fixed true stellar mass) applied. Note that deconvolution applied to noisy data (e.g., SMFs) amplifies the noise, 
which is why features in the left panel appear smoother than those in the right panel. 



TABLE 7 

Empirical Fit to Cosmic Star Formation Rates 



Publication 


zo 


A 


B 


C 


This Work 


1.243 


-0.997 


0.241 


0.180 


Hopkins & Beacom (2006) 


0.840 


-1.311 


0.085 


0.143 



NOTE. — Fits according to Eq. lFll for the collected cosmic star formation rate data in Table [4] ("This Work") and the data in Hopkins & Beacom (2006), all 
corrected to a Chabrier (2003) IMF. The quality of the fits may be seen in Fig. [2] The constant C has units of M© yr~' comoving Mpc~ 3 ; all other constants are 
unitless. Systematic errors are estimated in Table|3] 



within 0.02 dex in the vast majority of cases. Hence, to avoid overfitting bias, we reduce the magnitude of errors by 0.02 dex at 
each data point before calculating the \ 2 value used with our MCMC method. At z < 0.2, this has the impact of decoupling the 
systematic uncertainty parameters from the detailed shape of the SMF; at higher redshifts, there is no effect, as the typical size 
scale of the errors for observational data is already 0. 1-0.2 dex. 

THE REDSHIFT EVOLUTION OF THE STELLAR MASS - HALO MASS RELATION FROM ABUNDANCE MATCHING 

As a check of the accuracy of our analysis pipeline, it is instructive to compare with a direct abundance matching approach. In 
Fig- EU we show results for both direct abundance matching (left panel) and abundance matching including an assumed intrinsic 
(£) and measurement (a) scatter that match our assumed fiducial defaults. Notably, several features present in our main analysis 
are evident in this simpler approach: 

1. The SMHM relation for low-mass halos (M/, < 10 n Mo) at z = deviates from a power law. 

2. The halo mass at which the integrated star formation is most efficient rises by 0.5 dex from 10 12 M Q at z = to 10 12 5 M at 
z = 2. 

3. Above z = 2 it is difficult to locate the halo mass at which the peak integrated efficiency occurs because surveys at high 
redshifts do not have sufficient volume to probe the exponential tail of the stellar mass function (if one exists). If the UV 
luminosity function is any indication, however, the z = 8 results suggest that the peak may occur at lower masses above 
z = 2. 

4. There appears to be a turnaround in the SMHM relation near z ~ 2. The integrated efficiency of M/, ~ 10 11 5 M Q halos falls 
from z = to z = 2, but then increases at higher redshifts. 

At high redshifts, it is difficult to tell what, if any, evolution there is in the shape of the SMHM relation. However, it is clear that 
either the efficiency or the characteristic halo mass scale is evolving. If we constrain the shape parameters in our fitting formula 
for the SMHM relation (a, 8, and 7) to the z = values, we find that high redshifts strongly prefer a halo mass scale that scales as 
a power law with redshift, but have comparatively little implied evolution in the efficiency at high redshifts. For that reason, when 
we extrapolate to high redshifts unconstrained by data in our main analysis (z > 8.5), we assume that e, a, S, and 7 do not evolve 
beyond z = 8.5, but that the characteristic halo mass M; continues to evolve according to the redshift scaling we have adopted. 




10 



10 



FIG. 22. — Plots of the s pecifi c star formation rate data assembled in this study (see Table[6) for galaxies of stellar mass 10 9 5 M Q (top) and 10 ia5 M Q (bottom). 
Double power-law fits (Eq. lFU to th e data sets are taken from Tablets] Error bars are not shown since authors often dramatically underestimate the true magnitude 
of their systematic errors (see j|3.3t . 



TABLE 8 

Empirical Fits for Specific Star Formation Rates 



Stellar Mass 


zo 


A 


B 


C 




10 s Af 


1.018* 


-1.130 


-0.071 


3.659 x 10 


-9 


1O 95 M 


1.018* 


-1.259 


-0.053 


2.001 x 10 


-9 


10 l0 M o 


1.018 


-1.307 


-0.085 


1.462 x 10 


-9 


10 la5 M Q 


1.096 


-1.324 


-0.118 


9.032 x 10" 


10 



NOTE. — Fits a ccord ing to Eq. lFll for the collected specific star formation rate data in Table|6] all corrected to a Chabrier 12003) IMF. The quality of the fits 
may be seen in Fig. 1221 Asterisks denote that the parameter in question (zo) was unconstrained by the data at the given stellar mass and so was fixed to the value 
of the parameter at 10 i0 Mq . The constant C has units of yr ; all other constants are unitless. Systematic errors are estimated in Table|9] 



TABLE 9 

Error Estimates for Specific Star Formation Rates 



Redshift Range 


1O 9 M 


10 9i M o 


10 1U M Q 


10 lui M Q 


0.025-0.5 




0.39 


0.31 


0.29 


0.5-0.9 




0.32 


0.31 


0.27 


0.9-1.7 




0.30 


0.19 


0.29 


1.7-3 






0.24 


0.35 


3-8 


0.22 


0.22 


0.21 


0.31 



Average 0.28 



NOTE. — The standard devi ation in dex between published data p oints (from references in Table[6) in redshift bins after correction for the evolution of the 
specific star formation rate (Eq. lFU . For reasons mentioned in j|3.3l the average does not include the estimates in the z = 3-8 redshift range. The error in the 
error estimate for each bin is of order 45%; the error on the average error is of order 12%. Errors are not estimated when there are fewer than four references in a 
given stellar mass and redshift range. 

CALCULATING UNCERTAINTIES IN COSMIC AND SPECIFIC STAR FORMATION RATES 

For cosmic star formation rates (SFRs), we calculate the inter-publication variance in five redshift bins with 16 data points 
each. To take into account variation of the cosmic SFR in bins, we subtract a double-power law fit to the data points in Fig. [2J 

CSFR ^ = W^W^y (F1) 

where A, B, C, and zo are constants given in Table [7j adding more parameters did not make the quality of the fit substantially 
better. Our results are given in Table |5j we find the average systematic error to range from 0.13 dex at z = to 0.27 dex for z > 3 . 

For SSFRs (Tabled, we interpolate data for each publication's tabulated redshift ranges to obtain (if possible) estimates of the 
SSFR at four stellar masses (10 9 , 10 95 , 10 10 and 10 10 - 5 M Q ). We divide the SSFR data into the same five redshift bins as for the 
cosmic SFR data and estimate the systematic errors from the inter-publication variance. To exclude the effects of evolution in 
the SSFR over a redshift bin, we again subtract a double-power-law fit (Eq. lFll) with parame ters determined separate ly for each 
stellar mass in Table [HI the quality of the fit is shown in Fig. [22] It is worth remarking that I Weinmann et aT] d201 ll) conduct a 
similar analysis of SSFRs at M* = 10 9 5 M Q and find a "plateau" at high redshifts; while the data are consistent with a plateau, 
they are also consistent with not having a plateau: indeed, any range of shallow slopes at high redshifts would be consistent with 
the data collected here. 

Our results for the inter-publication variance as a function of stellar mass and redshift are shown in Table [9] Because only 
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FIG. 23. — Left: A comparison between the Tinker et al. ( 2008) mass function and two simulations with different initial conditions (ZA: Zel'dovich Approx- 
imation; 2LPT: second-order Lagrangian Perturbation Theory). Specifically, the plot shows the ratio of the cumulative central halo mass function to 10" 5 Mq 
(<E>(10" 5 Mq)) in the simulations to that for the Tinker MF. Small mass-dependent corrections (Eq. lGU were made to renormalize the Consuelo mass function to 
match the Tinker MF at z = 0, but no such corrections were necessary for Bolshoi. The Tinker MF, which was calibrated using ZA simulations, agree s remarkably 
well with Bolshoi o ver a wide range of redshifts; however, it is noticeably different from the Consuelo simulation for z > 3. As several authors ( Cr occe et al.1 
200 6; Je nkins 2010) have shown improved convergence for 2LPT over ZA initial conditions, we use the 2LPT normalization as the basis for our corrective fit 
(Eq. lG2t . Right: the mass-dependence of the correction to the Tinker cumulative mass function from the Consuelo simulation. Solid lines show direct results 
from the simulation; dashed lines show the fit in Eq. lG3l The dot-dashed line shows the fit in Eq . IG3 l extrapolated down to <3?(10 10 Mq). 

UV surveys can probe low stellar masses at high redshifts (z > 3), the published SSFRs for low stellar masses at high redshift 
are more in agreement (0.22 dex scatter) than for low stellar masses at low redshift (0.32 dex scatter)! This strongly suggests 
that the inter-publication variance in these bins does not represent the true systematic uncertainties. Excluding high-redshift bins 
(z > 3), there are no obvious systematic trends in the inter-publication variance, either with redshift or stellar mass. Indeed, a 
one-way analysis of variance (ANOVA) test delivers no strong evidence that the inter-publication variances differ across bins 
(p = 0.6). This conclusion is not affected by merging bins by stellar mass or by redshift, and the size of the effect is too small 
to be measurable without adding more data. We adopt the fiducial model of constant errors across stellar masses and redshifts, 
where the error is determined by merging all redshift and stellar mass bins except for the high-redshift bins (z > 3) which are 
known underestimates, we would estimate the average systematic error in SSFRs to be 0.28 dex. 

Finally, we mention a commonly-overlooked systematic error in optical SFRs, which comes from the bias introduced by 
separately fitting the dust model for each galaxy. As the dust model for each galaxy can never be derived exactly, individual 
galaxy SFRs will have some scatter relative to the true galaxy SFRs. One way to estima te the shape of this scatter is to compare 
two different SFR estimates with different assumptions; for example, Sa lim et all (12007b finds log-normal scatter consistent with 
intrinsic scatter of 0.3 dex. In calibrating one estimator to another (or to observations), authors typically make sure that the median 
offsets are zero; however, they often fail to notice that the expected mean offset is nonzero because the log-normal distribution is 
skewed in linear space. Specifically, a log-normal distribution with median value 1 and scatter a (in dex) has an expected mean 
value of 



(x)=exp(0.5(oinl0r). 



(F2) 



For a log-normal distribution with 0.3 dex scatter, this corresponds to an overall bias of +0.1 dex. This correction does not apply 
to UV studies that apply a median dust correction to all galaxies or to those that already correct for this bias. 

CORRECTIONS TO AND UNCERTAINTIES IN THE TINKER MASS FUNCTION 



The halo mass function for distinct halos presented by lTinkefet al. (2008) was calibrated for a wide array of simulations from 
z = to z = 2.5. The present work requires a mass fu nction from z = to z > 8 including satellite halos, so some additional 
corrections must be made. In addition, as new studies dReddick et al.H2012T) suggest that abundance matching to halo properties 
at their peak value before stripping offers a better match to z = clustering data, calibration for the quantities of interest using 
A^peak must also be applied. 

To test the high-redshift calibration of the mass function in iTinker' et al. (2008), we integrate the mass function in two simu- 
lations (Bolshoi and Consuelo) to 10 1L5 M Q and compare to the Tinker mass function for the appropriate cosmology. A mass- 
dependent correction is necessary to normalize the z = mass function for Consuelo to the Tinker mass function to account for 
slight incompleteness in the Consuelo mass function at this level: 



log 



10 



^Consuelo 



(M) 



AM) 



= -0.0008-0.042 



M 



10 1L51 M ( 







-0.69 



(Gl) 



The net effect of this correction on $(10 11 M Q ) is an increase of 0.03 dex (7%) compared to the raw s imulation result for 
Consuelo. No correction is necessary for the Bolshoi simulation. As shown in Fig. [23] the redshift fit given in (Tinker et al.l d2008) 
matches Bolshoi to within 0.02 dex up to a redshift of 5-6. However, Consuelo diverges from the Tinker mass function by 0.06 
dex (15%) already at z = 3. This discrepancy is due to the different initial conditions calculations used in the two simulations: 
Consuelo uses second-order Lagrangian Perturbation Theory (2LPT), whereas Bolshoi and the Tinker MF use the Zel'dovich 
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FIG. 24. — Best fits to the evolution of the sat ellite fractions (Eqs. lG7IG10t in Bolshoi for satellite accretion masses (left panel) and satellite peak masses (right 
panel), according to the simplified model of Eq. lGol The colored lines correspond to the fitted satellite fractions; the black dashed lines correspond to the satellite 
fractions in Bolshoi. 



approximation (1LPT ). The Zerdovich approximation has been previously found to underestimate the early nonlinear collapse 
of overdense regions dCrocce et al.ll2006tlJenkinsll2010l) . which is reconfirmed by the comparison in Fig. [23] For that reason, we 
prefer to use the correction from Consuelo, parametrized as follows: 

/ tfwqO^Me) \ _ 0.144 
810 V$Tinker(lO 11 - 5 M ) y ' 1 + exp [14.79(a- 0.213)] ' 

As shown in Fig. [23] this fit matches the correction from the Consuelo simulation to within 1%. However, as also shown in Fig. 
|23l there is a slight mass dependence to the correction; the magnitude of the correction is larger for more massive halos, as their 
nonlinear collapse begins earlier. We parametrize the mass dependence as follows: 



f $true(M) \ / ^(IO'^Mq) \ / M \ ^ 

l0gl ° { W^) J = 10810 W^M~ Q ) (G3) 



This function is well-behaved, even when extrapolated down to low masses, as shown in Fig. [23] 

Parametrizing the corrections as we have done introduces several uncertainties. The use of only a single simulation to calibrate 
the mass function introduces sample variance uncertainties, but these are in fact small for a box as large as Consuelo (420 Mpc 
/r 1 ); at low redshifts, the median offset from the Tinker ma ss function is 0.2% (see Eq. lGU . and even at redshift 6, the Poisson 
and sample variance errors combined are on the order of 1% dTrenti & Stiavellil2008l) for $(10 11,5 M Q ). Much larger uncertainties 
come from the use of a mass-dependent incompleteness correction to the Consuelo mass funct ion, because uncertainties in this 
correction are degenerate with the mass dependence of the normalization correction in Eq. IG3I If the incompleteness correction 
needed to be larger at high redshift, the evidence for the mass-dependence in the normalization correction would weaken. We can 
take this as an estimate of the additional systematic uncertainties on the high-redshift mass function. In terms of the impact on the 
halos in which high-redshift galaxies (z = 7-8) are expected to reside, the number density of such objects would be underestimated 
by at most 0.05 dex (12%), according to Fig. [23] if there were no mass-dependence in the normalization correction. This is tiny 
in comparison to the uncertainties in the stellar masses and star formation rates at high redshift. Indeed, even if the real universe 
corresponded exactly to the Tinker mass function, the total error introduced by using the corrections in this section (0. 1 dex, or 
26%) would again be dwarfed by uncertainties in measuring galaxy stellar content. 

In terms of the satellite fraction, we determine fits using a simple, physically-motivated model. It is an oft-observed feature 
of dark matter substructure that the number of satellite halos at a given mass scales approximately with the mass of the host — 
the number density of satellites in hosts of mass M c is proportional to M c <j> c (M c ), where <fi c (M) is the mass function for central 
halos. Another feature of dark matter substructure is that the number of satellites at a given satellite-to-host mass ratio is largely 
self-similar across host halo masses; this requires that the proportionality must scale inversely with the satellite mass. Therefore, 
we can write down a simple formula for the satellite mass function, (j> s : 



C(a) f°° 

MM) ~ 4r / M c (/> c (M c )d\og 10 M c (G4) 

IM 



M 



where the proportionality factor, C(a), is a function of scale factor only. The low-mass limit of this function takes on an especially 
appealing form; as <j> c (M c ) <~ <fioM~ l foiM c is less than the exponential cutoff scale, M cuto ff, we find: 

MM) ^C'"Vo^log 10 M t . 
0cC*O ^oM" 1 
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which can be simplified to 



'C(a)log 10 — n — (G6) 



C (M) v elu V M 



Thus, the satellite fraction gro ws o nly logarithmically in the low-mass limit. 

However, the picture in Eq. |G4| is not entirely complete; because they have later assembly times, massive clusters are more 
likely to be undergoing major mergers than smaller halos. However, these major mergers also become absorbed more quickly, as 
the effects of dynamical friction are more pronounce d in massive halos. These two (opposite) effects combine to give a satellite 
fraction slightly larger at the high-mass end than Eq. IG4I woul d pre dict; moreover, this excess evolves nontrivially with redshift. 
For that reason, we use the slightly more flexible form of Eq. IG6I and allow great freedom in the parametrization for both C(a) 
andM cutoff (a). 

Fitting the satellite fraction to a simulation is made especially tricky by the fact that simulations are incomplete for satellites 
at higher masses than for central halos, even before selecting satellites on M pea k or M acc . For that reason, Consuelo is unsuitable 
for determining the fraction of all but the most massive satellites (M ~ 10 13 M Q ), which provides an especially poor handle on 
constraining the redshift evolution of the fit. Bolshoi fares better, being complete for satellites down to 10 n M Q ; yet, given its use 
of ZA initial conditions, it may not seem especially suitable for high-redshift calibration, either! 

Yet, there are two arguments in its favor. First, the satellite fraction is low at high redshift, so even large relative errors in the 
satellite fraction will translate into small errors on the overall mass fraction. Second, the influence of the initial conditions on the 
satellite fraction are largely through the effects of assembly bias; however, recent work has shown that satellite richness (when 
tallied above a fixed cut in accretion or peak mass) is almost uncorrected with formation time (Wu et al., in prep.). The impact 
of the initial conditions on the satellite fraction is thus minor compared to the effect of the mass function normalization. 

Therefore, using the Bolshoi simulation, we fit the satellite fraction from z = to z = 10; when satellite masses are determined 
by M acc , we find: 

Zog 10 (C(a)) = -2.69+11.68a-28.88a 2 +29.33a 3 - 10.56a 4 (G7) 
log 10 (M cutoff (a)) = 1 1 .34 + 8 34a - 0.36a 2 - 5 .08a 3 + 0.75a 4 (G8) 

and find the following parametrization for the satellite fraction when satellite masses are determined by M pea k: 



log 10 (C(a)) = -1.91+6.23a-15.07a 2 +15.02a 3 -5.29a 4 (G9) 
log 10 (M cutoff (a))= 10.66+ 15.93a-21.39a 2 +18.20a 3 -8. 21a 4 (G10) 

Both of these fits in comparison to sat ellite fractions in B olshoi are shown in Fig.[24l they are accurate to within the estimated 
uncertainties in the mass function r5%. lTinker eTai]l2008h . 

HALO MASS ACCRETION RATES 

The most commonly used fitting formula for halo mass accretion rates is that proposed by Wec hsler et al.l d2002l) . which 
suggests that halo mass growth is an exponential function in redshift (M(z) oc exp(-cz)). As simulations improved, deviations 
fr om this simple f ormul a were found, which led authors to propose more complicated formulas; e.g., M(z) oc (1 +z) 6 exp(-cz) 
in McBrid eet alj d2009l) . However, little attention has been given to the high-redshift behavior of these formulas, with the 
predictable result that they fail to accurately capture median mass accretion histories for halos as recently as z ~ 3 (see Fig. [25] 
left panel)Q 

This limitation a problem for this study. To motivate our choice of improved fitting formulas, we note an issue with all proposed 
fitting formulas to date, which is their exponential dependence at high redshifts. That is to say, the mass ratio of two halos Mi 
and M% will asymptote to 

M,(z) M,(0) 
M 2 (z) M 2 (0) 

This is problematic because more massive halos usually have steeper exponential dependencies (i.e., Mi < M2 implies c\ < C2). 
If M\ (z = 0) < M2(z = 0), these formulas then predict that eventually the progenitor histories will cross and Mi (z) > M2(z). Using 
tfae lMcBride et al.l (|2009) formula, for example, with the MultiDark simulation, we find that it predicts the progenitors of 10 15 M Q 
halos to be less massive than the progenitors of 10 14 75 M© halos as early as z = 5. However, this is not the case: in simulations, 
the hierarchy of average or median progenitor masses is maintained, even though the distance between them may decrease (see, 
for example, the right panel of Fig. |25j. 

We therefore parametrize the mass accretion history of a single halo mass, and then parametrize the ratios of progenitor 
masses for other halo masses, which entirely avoids the problem of crossing progenitor mass histories. The mass bin we choose 
(10 13 15 - 10 13 40 Mq) is well-resolved in all three dark matter simulations at our disposal, and sufficient statistics are available 

4 Note that, for individual halos, as opposed to ensemble medians or aver- ... c . ,. ., . . 

, . , ' . . . ,.,,.„, .. , ^ r . , , , „. , . , cision is not as necessary tor cases when individual stochasticity in merger 

ages, the| McBndeetalj formula provides a reasonable fit: high pre- r ^ fa mgher deviations than those shown m the bottom . left panel 

ofFig.f25] 
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FIG . 25 . — Fits to median mass accretion histori es (M v j r ) for central ha los in Bolshoi, MultiDark, and Consuel o simulations. Commo n functional forms for the 
mass accretion rates, such as the exponential fit of Wechsler et al. (2002) or the exponential plus power-law fit of McBride et al. 12009) fail to accurately capture 
mass accretion histories beyond z = 2. Moreover, at high redshifts, these common functional forms imply unphysical behavior — i.e., that progenitors of 10 15 Mq 
halos are on average less massive than progenitors of 10 13 Mq halos (see text), which is not at all the case. The left panel shows mass accretion histories for 
jq13.15 to if) 13 40 Mq halos in Bolshoi and Consuelo, along with residuals (linear scale). Bolshoi and Consuelo have different cosmologies, so it is expected that 
the mass accretion rates of halos differ slightly. The right panel shows progenitor mass ratios (compared to progenitors of 10 13 Mq halos) for halos of mass 10" 
to 10 15 Mq at z = for the Bolshoi, MultiDark, and Consuelo simulations, along with residuals. In all cases, our fits to the mass accretion histories of halos for 
the Bolshoi cosmology are accurate to 5% on average; the only exception is for the highest mass bins in the Bolshoi and Consuelo simulations, which contains 
only a few halos. 



even in the smallest simulation (Bolshoi) to be able to reduce scatter to an acceptable level. Our resulting fit is: 

M(Mo,z)=M B (z)10 /(Mo ' z) (H2) 

M 13 (z)=10 13276 (l+z) 3 00 (l + |r 611 exp(-0.503z)M Q (H3) 

Wo,z)=1oJ^)J^ (H4) 



M 13 (0); g(Mo,^) 

g(M ,a) = 1 + exp(-4.65 1 (a - a (M )) (H5) 



flo(M ) = 0.205 -log, 



10 9649 M, 







M 



0.18 

+ 1 



(H6) 



where M(Mq,z) gives the median virial mass dBryan & Nor man 1998) for progenitors of halos with mass Mq at z = 0. Despite its 
complexity, this seven-parameter fit (three parameters for M\^{z), one parameter for g{M(,,a), and three parameters for a^Mo)) 
provides an excellent fit to mass accretion histories for a wide range of halo masses and redshifts, as shown in Fig. [25] 

SUBHALO MERGER/DISRUPTION RATES 

Halo mergers are important for our study because they impact the rate at which stellar mass is transfered from low-mass galaxies 
to high-mass galaxies as well as the rate at which stars are ej ected into the intraclu ster light (ICL). While there have been many 
studies of the rate at which halos become subhalos (see, e.g.. lFakhouri et afll2010l) in dark matter simulations, these studies are 
less relevant to our work because it is expected that sat ellite galaxies should last at least as long as the corresponding subhalo that 
hosts them (e.g. lYang et al.ll20l2l iReddick et alj|2012l) . Studying the subhalo destructio n rate is more difficult because of issues 
with halo finding completeness when subhalos pass close to the cent ers of their hosts (|Behroozi et al1l201 1 alibi) , and relatively 
few authors have attempted to correct for this issue in their work (c.f., Wetzel & Whitell2010F ~ 

However, with the advantage of a phase-space halo finder and a merger tree algorithm that corrects for such inconsistencies, 
we are in a unique position to provide a calibration of the merger rate (or disruption rate) of subhalos into host h alos all the way 
back to the earliest redshifts in simulations. We find that the functional form proposed in lFakhouri et al.1 (120101) is an excellent 
fit to the subhalo disruption rate; however, the parameter values and redshift evolution are different due to the different merger 
definition, cosmology, and mass definition used. 

Letting M/, be the mass of the host halo and 6 be the ratio of subhalo peak mass to host mass, we find that the disruption rate 
(expressed in terms of the decrement in number of subhalos per unit host halo per unit redshift per log interval in subhalo mass 
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FIG. 26. — The rate of subhalo disruption in host halos in the Bolshoi simulation with mass and redshift dependencies (X(M, z)) scaled out. Results are presented 
in terms of the number of subhalos disrupted per host halo per unit redshift per log interval in satellite mass ratio. The subhalo merger mass ratio (9) in these 
plots is the ratio of the subhalo peak mass (i.e., the highest mass any progenitor of the subhalo ever had) to the host halo mass. The left panel shows the scaled 
subhalo disruption rate at z ~ over a range of host halo masses, demonstrating excellent convergence. The red dashed line in the left panel represents the 
average of scaled subhalo disruption rates as a function of subhalo mass ratio 6, excluding halos for which subhalos of the given mass ratio would be below the 
mass resolution limit of the simulation (1O'°M0). The right panel shows the mass-averaged scaled subhalo disruption rate as a function of redshift, as well as the 
fit we adopt in this work. The fit shows excellent conformance to the simulation data, except for very low subhalo mass ratios, where stochasticity in evaluating 
halo masses at low particle numbers artificially inflates peak halo masses (see text). 
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We show the quality of the fits for the peak subhalo mass in Fig. [26] For ease of comparing simulation data to the fit, we have 
scaled out host mass and redshift dependencies according to Eqs. llTllBl so that all plotted lines are directly comparable to the 
subhalo disruption rate forM = 10 I2 M Q at z = 0. 

Fig.[26lshows that scaling in both mass and redshift of the subhalo disruption rate has been successfully captured in Eqs. lTITlBl 
The remaining deviations are due to Poisson noise in the disruption rates for halo masses where Bolshoi has few halos (especially 
in the highest mass bin) and stochasticity in halo mass estimations for low-mass halos. Specifically, as the peak mass captures the 
highest mass ever recorded in a halo's mass accretion history, greater variance in the halo mass estimator translates to artificially 
increased peak masses for low-mass halos. This manifests as a slight upturn in the subhalo disruption rate close to the halo mass 
resolution limit of the simulation. 



FITS TO STAR FORMATION RATES AND HISTORIES UNDER DIFFERENT MODELING CONSTRAINTS 

We consider four models with alternate sets of assumptions to the main analysis in this paper: 

1 . A model ("No Systematics") that excludes all nuisance parameters related to errors in converting luminosities into stellar 
masses (i.e., [i = n = 0, Tfcai = oo) and to incompleteness (A = 0). 

2. A mode l ("No z > 6 C onstraint s") that excludes all da ta constr aints at high redshift (z > 6), including data in lBouwens et alJ 
(1201 11): iBradlev etalJ (120121) : ISchaerer & de Barrosl (120101) : iMcLure et all (120111) . and the data points above z = 6 in 
iBouwens et al.l (l2012h . 

3. A model ("No CSFR Constraints") that excludes all constraints from observed cosmic star formation rates. 

4. A model ("No SFR Constraints") that excludes all constraints from observed star formation rates (both CSFRs and SSFRs) 
as well as excluding all systematic nuisance parameters as in the first alternate model ("No Systematics"). 
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FIG. 27. — Comparisons to observational data for the best fits of alternate models. Top left: comparisons to observed stellar mass functions for the model 
without systematic uncertaint ies (i.e., fj, = 0, k — 0, T) ca ) = oo, no incomp l eteness). Top right: c o mparisons to stellar mass fun c tions for a model exc luding all 
constraints at z > 6 (i.e., frorr[Bouwens et al. 1 201 1); Bradley et al. (2012); Bo uwens etal] 12012): Sch aerer & de Barrosl )2010l) ; IMcLure et alj 1201 II) ). Middle 
left: comparisons to observations of the cosmic star formation rate for the two models already mentioned, along with a model that excludes all cosmic star 
formation rate constraints ("No CSFR"; using SSFRs only), and a model that excludes all star formation rate constraints and systematic effects ("No SFR"). 
Middle-right and bottom panels: comparisons with observed SSFRs for all five models. 

Comparisons with observed data are presented in Fig. 127] For ease of comparison across mass ranges, we show comparisons 
between derived baryon conversion efficiencies and stellar mass histories for the main model in this paper and the alternate 
models in Figs.[28land|29l 

First, it may be verified that the model including no systematics succeeds remarkably well in matching all three kinds of 
observations; while it is potential ly so mewhat high at low redshifts for the CSFR, this is not unexpected on account of issues 
with the galaxy / ICL definition Q2.41 >. The No Systematics model has much tighter error bars on the stellar mass to halo mass 
ratios (not shown) than the full model at all redshifts, which suggests that systematics remain the single largest source of errors 
in determining stellar masses. Fig.|28]is comparable to Fig.|5]for the No Systematics model, except that at high redshifts and low 
halo masses, there is an upturn in the SFR that is not present in the full model. 

The model excluding z > 6 data gives almost identical results as the full model for z < 6, which is encouraging. The "predic- 
tions" for high-redshift stellar mass functions are somewhat lower than in the best-fitting model, which results in a higher SFR 
(i.e., rate of stellar mass growth) over the period from z ~ 7 to z ~ 6. 

The model excluding all cosmic star formation constraints is remarkably similar to the full model in our analysis. This is as 
expected: the specific star formation rate in combination with the stellar mass function gives a constraint on the total amoun t of 
star formation. Yet, the CSFR data has tighter error bars than the SSFR in combination with the SMF, at least for z < 1 (see jj |3.2| 
and l3.3l >. so it provides slightly better constraints at those redshifts. 

Excluding all star formation rate data results in a good fit to the CSFR at high redshifts, where growth in stellar mass is rapid — 
i.e., where stellar mass functions at successive redshifts are significantly different even despite large errors. At lower redshifts, 
the predicted cosmic star formation rate for the best-fitting model is discrepant from observations. At these redshifts, the growth 
in stellar mass is poorly constrained because the intrinsic stellar mass functions are changing comparatively little, especially from 
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FIG. 28. — Derived baryon conversion efficiencies and stellar mass histories for three alternate models (excluding all systematic uncertainties, excluding 
constraints at z > 6, and excluding constraints from cosmic star formation rates); the shaded bands show the one-sigma distributions from the main model in this 
paper, and the solid lines show results from the alternate models. 
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FIG. 29. — Derived baryon conversion efficiencies and stellar mass histories for one alternate model (excluding all systematic uncertainties and constraints from 
star formation rates); the shaded bands show the one-sigma distributions from the main model in this paper, and the solid lines show results from the alternate 
model. 

Z = 0.5 to z = 0. A small difference in the choice of stellar mass function evolution at low redshifts therefore has a much larger 
impact on the relative error in star formation rates. Thus, although the stellar mass to halo mass ratios are similar to the No 
Systematics model (and the abundance matching model in Appendix |EJ, the star formation rates derived using this model are 
unreliable. 



